From 8176291618baac5e7d6e49e367687eac051d98d6 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 24 Nov 2015 17:46:31 +0100 Subject: [PATCH 01/31] Maybe repaired perturb_buffer_by_mono --- plugins/Perturbation/perturbation.template.f | 49 ++++++++++++++------ 1 file changed, 34 insertions(+), 15 deletions(-) diff --git a/plugins/Perturbation/perturbation.template.f b/plugins/Perturbation/perturbation.template.f index fa01cc99..b41c7685 100644 --- a/plugins/Perturbation/perturbation.template.f +++ b/plugins/Perturbation/perturbation.template.f @@ -54,8 +54,6 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c cycle endif - integer :: degree - call get_excitation_degree(HF_bitmask,buffer(1,1,i),degree,N_int) call pt2_$PERT(psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, & c_pert,e_2_pert,H_pert_diag,Nint,N_minilist,n_st,minilist,idx_minilist,N_minilist) @@ -72,7 +70,8 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c end -subroutine perturb_buffer_by_mono_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint) + +subroutine perturb_buffer_by_mono_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint,key_mask,fock_diag_tmp) implicit none BEGIN_DOC ! Applly pertubration ``$PERT`` to the buffer of determinants generated in the H_apply @@ -81,20 +80,41 @@ subroutine perturb_buffer_by_mono_$PERT(i_generator,buffer,buffer_size,e_2_pert_ integer, intent(in) :: Nint, N_st, buffer_size, i_generator integer(bit_kind), intent(in) :: buffer(Nint,2,buffer_size) + integer(bit_kind),intent(in) :: key_mask(Nint,2) + double precision, intent(in) :: fock_diag_tmp(2,0:mo_tot_num) double precision, intent(inout) :: sum_norm_pert(N_st),sum_e_2_pert(N_st) double precision, intent(inout) :: coef_pert_buffer(N_st,buffer_size),e_2_pert_buffer(N_st,buffer_size),sum_H_pert_diag(N_st) double precision :: c_pert(N_st), e_2_pert(N_st), H_pert_diag(N_st) - integer :: i,k, c_ref + integer :: i,k, c_ref, ni, ex integer, external :: connected_to_ref_by_mono logical, external :: is_in_wavefunction + integer(bit_kind) :: minilist(Nint,2,N_det_selectors) + integer :: idx_minilist(N_det_selectors), N_minilist + + integer(bit_kind) :: minilist_gen(Nint,2,N_det_generators) + integer :: N_minilist_gen + logical :: fullMatch + logical, external :: is_connected_to + + + ASSERT (Nint > 0) ASSERT (Nint == N_int) ASSERT (buffer_size >= 0) ASSERT (minval(sum_norm_pert) >= 0.d0) ASSERT (N_st > 0) - do i = 1,buffer_size - + + call create_minilist(key_mask, psi_selectors, miniList, idx_miniList, N_det_selectors, N_minilist, Nint) + call create_minilist_find_previous(key_mask, psi_det_generators, miniList_gen, i_generator-1, N_minilist_gen, fullMatch, Nint) + + if(fullMatch) then + return + end if + + + do i=1,buffer_size + c_ref = connected_to_ref_by_mono(buffer(1,1,i),psi_det_generators,Nint,i_generator,N_det) if (c_ref /= 0) then @@ -105,19 +125,18 @@ subroutine perturb_buffer_by_mono_$PERT(i_generator,buffer,buffer_size,e_2_pert_ cycle endif - integer :: degree - call pt2_$PERT(buffer(1,1,i), & - c_pert,e_2_pert,H_pert_diag,Nint,N_det_selectors,n_st) + call pt2_$PERT(psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, & + c_pert,e_2_pert,H_pert_diag,Nint,N_minilist,n_st,minilist,idx_minilist,N_minilist) do k = 1,N_st - e_2_pert_buffer(k,i) = e_2_pert(k) - coef_pert_buffer(k,i) = c_pert(k) - sum_norm_pert(k) += c_pert(k) * c_pert(k) - sum_e_2_pert(k) += e_2_pert(k) - sum_H_pert_diag(k) += H_pert_diag(k) + e_2_pert_buffer(k,i) = e_2_pert(k) + coef_pert_buffer(k,i) = c_pert(k) + sum_norm_pert(k) = sum_norm_pert(k) + c_pert(k) * c_pert(k) + sum_e_2_pert(k) = sum_e_2_pert(k) + e_2_pert(k) + sum_H_pert_diag(k) = sum_H_pert_diag(k) + H_pert_diag(k) enddo - enddo + enddo end From d05c851ba62401d4ebee57a2f2ce7a399e56e6e7 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 25 Nov 2015 01:08:31 +0100 Subject: [PATCH 02/31] cleaning --- src/Integrals_Bielec/map_integrals.irp.f | 22 +++--- src/Integrals_Bielec/mo_bi_integrals.irp.f | 3 +- src/Utils/map_module.f90 | 91 ---------------------- 3 files changed, 13 insertions(+), 103 deletions(-) diff --git a/src/Integrals_Bielec/map_integrals.irp.f b/src/Integrals_Bielec/map_integrals.irp.f index 23deec02..f835f11f 100644 --- a/src/Integrals_Bielec/map_integrals.irp.f +++ b/src/Integrals_Bielec/map_integrals.irp.f @@ -318,6 +318,7 @@ double precision function get_mo_bielec_integral_schwartz(i,j,k,l,map) get_mo_bielec_integral_schwartz = dble(tmp) end + double precision function mo_bielec_integral(i,j,k,l) implicit none BEGIN_DOC @@ -356,36 +357,37 @@ subroutine get_mo_bielec_integrals(j,k,l,sze,out_val,map) call map_get_many(map, hash, tmp_val, sze) ! Conversion to double precision do i=1,sze - out_val(i) = tmp_val(i) + out_val(i) = dble(tmp_val(i)) enddo endif end -subroutine get_mo_bielec_integrals_existing_ik(j,l,sze,out_array,map) +subroutine get_mo_bielec_integrals_ij(k,l,sze,out_array,map) use map_module implicit none BEGIN_DOC ! Returns multiple integrals in the MO basis, all - ! i(1)j(1) 1/r12 k(2)l(2) - ! i for j,k,l fixed. + ! i(1)j(2) 1/r12 k(1)l(2) + ! i, j for k,l fixed. END_DOC - integer, intent(in) :: j,l, sze + integer, intent(in) :: k,l, sze logical, intent(out) :: out_array(sze,sze) type(map_type), intent(inout) :: map - integer :: i,k,kk,ll,m + integer :: i,j,kk,ll,m integer(key_kind),allocatable :: hash(:) integer ,allocatable :: pairs(:,:), iorder(:) PROVIDE mo_bielec_integrals_in_map allocate (hash(sze*sze), pairs(2,sze*sze),iorder(sze*sze)) kk=0 - do k=1,sze + out_array = 0.d0 + do j=1,sze do i=1,sze kk += 1 !DIR$ FORCEINLINE call bielec_integrals_index(i,j,k,l,hash(kk)) pairs(1,kk) = i - pairs(2,kk) = k + pairs(2,kk) = j iorder(kk) = kk enddo enddo @@ -404,8 +406,8 @@ subroutine get_mo_bielec_integrals_existing_ik(j,l,sze,out_array,map) do ll=1,kk m = iorder(ll) i=pairs(1,m) - k=pairs(2,m) - out_array(i,k) = (hash(ll) /= 0_8) + j=pairs(2,m) + out_array(i,j) = (hash(ll) /= 0_8) enddo deallocate(pairs,hash,iorder) diff --git a/src/Integrals_Bielec/mo_bi_integrals.irp.f b/src/Integrals_Bielec/mo_bi_integrals.irp.f index 98737e2a..1fa303b5 100644 --- a/src/Integrals_Bielec/mo_bi_integrals.irp.f +++ b/src/Integrals_Bielec/mo_bi_integrals.irp.f @@ -102,7 +102,7 @@ subroutine add_integrals_to_map(mask_ijkl) !$OMP mo_coef_transp, & !$OMP mo_coef_transp_is_built, list_ijkl, & !$OMP mo_coef_is_built, wall_1, abort_here, & - !$OMP mo_coef,mo_integrals_threshold,ao_integrals_map,mo_integrals_map,progress_bar,progress_value) + !$OMP mo_coef,mo_integrals_threshold,mo_integrals_map,progress_bar,progress_value) n_integrals = 0 allocate(bielec_tmp_3(mo_tot_num_align, n_j, n_k), & bielec_tmp_1(mo_tot_num_align), & @@ -315,7 +315,6 @@ IRP_ENDIF call ezfio_set_integrals_bielec_disk_access_mo_integrals("Read") endif - end diff --git a/src/Utils/map_module.f90 b/src/Utils/map_module.f90 index ecff478f..1c89355d 100644 --- a/src/Utils/map_module.f90 +++ b/src/Utils/map_module.f90 @@ -437,97 +437,6 @@ call omp_unset_lock(map%lock) end -subroutine map_update_verbose(map, key, value, sze, thr) - use map_module - implicit none - type (map_type), intent(inout) :: map - integer, intent(in) :: sze - integer(key_kind), intent(inout) :: key(sze) - real(integral_kind), intent(inout) :: value(sze) - real(integral_kind), intent(in) :: thr - - integer :: i - integer(map_size_kind) :: idx_cache, idx_cache_new - integer(cache_map_size_kind) :: idx - integer :: sze2 - integer(cache_key_kind) :: cache_key - integer(map_size_kind) :: n_elements_temp - type (cache_map_type) :: local_map - logical :: map_sorted -! do i = 1, sze -! print*,'value in map = ',value(i) -! enddo - - sze2 = sze - map_sorted = .True. - - n_elements_temp = 0_8 - n_elements_temp = n_elements_temp + 1_8 - do while (sze2>0) - i=1 - do while (i<=sze) - if (key(i) /= 0_8) then - idx_cache = ishft(key(i),map_shift) - if (omp_test_lock(map%map(idx_cache)%lock)) then - local_map%key => map%map(idx_cache)%key - local_map%value => map%map(idx_cache)%value - local_map%sorted = map%map(idx_cache)%sorted - local_map%map_size = map%map(idx_cache)%map_size - local_map%n_elements = map%map(idx_cache)%n_elements - do - !DIR$ FORCEINLINE - call search_key_big_interval(key(i),local_map%key, local_map%n_elements, idx, 1, local_map%n_elements) - if (idx > 0_8) then -! print*,'AHAAH' -! print*,'local_map%value(idx) = ',local_map%value(idx) - local_map%value(idx) = local_map%value(idx) + value(i) -! print*,'not a new value !' -! print*,'local_map%value(idx) = ',local_map%value(idx) - else - ! Assert that the map has a proper size - if (local_map%n_elements == local_map%map_size) then - call cache_map_unique(local_map) - call cache_map_reallocate(local_map, local_map%n_elements + local_map%n_elements) - call cache_map_shrink(local_map,thr) - endif - cache_key = iand(key(i),map_mask) - local_map%n_elements = local_map%n_elements + 1_8 - local_map%value(local_map%n_elements) = value(i) -! print*,'new value !' - local_map%key(local_map%n_elements) = cache_key - local_map%sorted = .False. - n_elements_temp = n_elements_temp + 1_8 - endif ! idx > 0 - key(i) = 0_8 - i = i+1 - sze2 = sze2-1 - if (i>sze) then - i=1 - endif - if ( (ishft(key(i),map_shift) /= idx_cache).or.(key(i)==0_8)) then - exit - endif - enddo - map%map(idx_cache)%key => local_map%key - map%map(idx_cache)%value => local_map%value - map%map(idx_cache)%sorted = local_map%sorted - map%map(idx_cache)%n_elements = local_map%n_elements - map%map(idx_cache)%map_size = local_map%map_size - map_sorted = map_sorted .and. local_map%sorted - call omp_unset_lock(map%map(idx_cache)%lock) - endif ! omp_test_lock - else - i=i+1 - endif ! key = 0 - enddo ! i -enddo ! sze2 > 0 -call omp_set_lock(map%lock) -map%n_elements = map%n_elements + n_elements_temp -map%sorted = map%sorted .and. map_sorted -call omp_unset_lock(map%lock) - -end - subroutine map_append(map, key, value, sze) use map_module implicit none From bef37a721c169111df36dcabe7fb9f9fd6c97c22 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 25 Nov 2015 02:01:40 +0100 Subject: [PATCH 03/31] Filter_connected accelerated --- src/Determinants/filter_connected.irp.f | 30 +++++++++---------------- 1 file changed, 11 insertions(+), 19 deletions(-) diff --git a/src/Determinants/filter_connected.irp.f b/src/Determinants/filter_connected.irp.f index 005aa020..1bf76dc4 100644 --- a/src/Determinants/filter_connected.irp.f +++ b/src/Determinants/filter_connected.irp.f @@ -130,9 +130,7 @@ subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx) do i=1,sze degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & popcnt(xor( key1(1,2,i), key2(1,2))) - if (degree_x2 > 4) then - cycle - else if(degree_x2 /= 0)then + if (degree_x2 <= 4) then idx(l) = i l = l+1 endif @@ -146,9 +144,7 @@ subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx) popcnt(xor( key1(2,1,i), key2(2,1))) + & popcnt(xor( key1(1,2,i), key2(1,2))) + & popcnt(xor( key1(2,2,i), key2(2,2))) - if (degree_x2 > 4) then - cycle - else if(degree_x2 /= 0)then + if (degree_x2 <= 4) then idx(l) = i l = l+1 endif @@ -164,9 +160,7 @@ subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx) popcnt(xor( key1(2,2,i), key2(2,2))) + & popcnt(xor( key1(3,1,i), key2(3,1))) + & popcnt(xor( key1(3,2,i), key2(3,2))) - if (degree_x2 > 4) then - cycle - else if(degree_x2 /= 0)then + if (degree_x2 <= 4) then idx(l) = i l = l+1 endif @@ -174,6 +168,8 @@ subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx) else + + integer, save :: icount(4) = (/0,0,0,0/) !DIR$ LOOP COUNT (1000) outer: do i=1,sze degree_x2 = 0 @@ -181,21 +177,17 @@ subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx) do m=1,Nint if ( key1(m,1,i) /= key2(m,1)) then degree_x2 = degree_x2+ popcnt(xor( key1(m,1,i), key2(m,1))) - if (degree_x2 > 4) then - cycle outer - endif endif if ( key1(m,2,i) /= key2(m,2)) then degree_x2 = degree_x2+ popcnt(xor( key1(m,2,i), key2(m,2))) - if (degree_x2 > 4) then - cycle outer - endif + endif + if (degree_x2 > 4) then + cycle outer endif enddo - if(degree_x2 /= 0)then - idx(l) = i - l = l+1 - endif + idx(l) = i + l = l+1 + icount(3) = icount(3) + 1_8 enddo outer endif From 1254176e9023a10c8cf9966e3a41ad060b7abc82 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 25 Nov 2015 03:02:19 +0100 Subject: [PATCH 04/31] Prefetching in map_integrals --- src/Integrals_Bielec/map_integrals.irp.f | 13 ++- src/Utils/map_module.f90 | 141 +++++++++++++++++++++-- 2 files changed, 140 insertions(+), 14 deletions(-) diff --git a/src/Integrals_Bielec/map_integrals.irp.f b/src/Integrals_Bielec/map_integrals.irp.f index f835f11f..3b737f5d 100644 --- a/src/Integrals_Bielec/map_integrals.irp.f +++ b/src/Integrals_Bielec/map_integrals.irp.f @@ -371,13 +371,16 @@ subroutine get_mo_bielec_integrals_ij(k,l,sze,out_array,map) ! i, j for k,l fixed. END_DOC integer, intent(in) :: k,l, sze - logical, intent(out) :: out_array(sze,sze) + double precision, intent(out) :: out_array(sze,sze) type(map_type), intent(inout) :: map integer :: i,j,kk,ll,m integer(key_kind),allocatable :: hash(:) integer ,allocatable :: pairs(:,:), iorder(:) + real(integral_kind), allocatable :: tmp_val(:) + PROVIDE mo_bielec_integrals_in_map - allocate (hash(sze*sze), pairs(2,sze*sze),iorder(sze*sze)) + allocate (hash(sze*sze), pairs(2,sze*sze),iorder(sze*sze), & + tmp_val(sze*sze)) kk=0 out_array = 0.d0 @@ -401,16 +404,16 @@ subroutine get_mo_bielec_integrals_ij(k,l,sze,out_array,map) call i2radix_sort(hash,iorder,kk,-1) endif - call map_exists_many(mo_integrals_map, hash, kk) + call map_get_many(mo_integrals_map, hash, tmp_val, kk) do ll=1,kk m = iorder(ll) i=pairs(1,m) j=pairs(2,m) - out_array(i,j) = (hash(ll) /= 0_8) + out_array(i,j) = tmp_val(ll) enddo - deallocate(pairs,hash,iorder) + deallocate(pairs,hash,iorder,tmp_val) end integer*8 function get_mo_map_size() diff --git a/src/Utils/map_module.f90 b/src/Utils/map_module.f90 index 1c89355d..24f5a328 100644 --- a/src/Utils/map_module.f90 +++ b/src/Utils/map_module.f90 @@ -496,13 +496,16 @@ subroutine cache_map_get_interval(map, key, value, ibegin, iend, idx) integer(cache_map_size_kind), intent(in) :: ibegin, iend real(integral_kind), intent(out) :: value integer(cache_map_size_kind), intent(inout) :: idx + double precision, pointer :: v(:) + integer :: i - call search_key_big_interval(key,map%key, map%n_elements, idx, ibegin, iend) - if (idx > 0) then - value = map%value(idx) - else - value = 0._integral_kind - endif +! call search_key_big_interval(key,map%key, map%n_elements, idx, ibegin, iend) + call search_key_value_big_interval(key, value, map%key, map%value, map%n_elements, idx, ibegin, iend) +! if (idx > 0) then +! value = v(idx) +! else +! value = 0._integral_kind +! endif end @@ -612,7 +615,7 @@ subroutine search_key_big_interval(key,X,sze,idx,ibegin_in,iend_in) istep = ishft(iend-ibegin,-1) idx = ibegin + istep - do while (istep > 32) + do while (istep > 16) idx = ibegin + istep if (cache_key < X(idx)) then iend = idx @@ -649,8 +652,8 @@ subroutine search_key_big_interval(key,X,sze,idx,ibegin_in,iend_in) endif enddo idx = ibegin - if (min(iend_in,sze) > ibegin+64) then - iend = ibegin+64 + if (min(iend_in,sze) > ibegin+16) then + iend = ibegin+16 !DIR$ VECTOR ALIGNED do while (cache_key > X(idx)) idx = idx+1 @@ -693,6 +696,126 @@ subroutine search_key_big_interval(key,X,sze,idx,ibegin_in,iend_in) end +subroutine search_key_value_big_interval(key,value,X,Y,sze,idx,ibegin_in,iend_in) + use map_module + implicit none + integer(cache_map_size_kind), intent(in) :: sze + integer(key_kind) , intent(in) :: key + real(integral_kind) , intent(out) :: value + integer(cache_key_kind) , intent(in) :: X(sze) + real(integral_kind) , intent(in) :: Y(sze) + integer(cache_map_size_kind), intent(in) :: ibegin_in, iend_in + integer(cache_map_size_kind), intent(out) :: idx + + integer(cache_map_size_kind) :: istep, ibegin, iend, i + integer(cache_key_kind) :: cache_key + + if (sze /= 0) then + continue + else + idx = -1 + value = 0.d0 + return + endif + cache_key = iand(key,map_mask) + ibegin = min(ibegin_in,sze) + iend = min(iend_in,sze) + if ((cache_key > X(ibegin)) .and. (cache_key < X(iend))) then + + istep = ishft(iend-ibegin,-1) + idx = ibegin + istep + do while (istep > 16) + idx = ibegin + istep + if (cache_key < X(idx)) then + iend = idx + istep = ishft(idx-ibegin,-1) + idx = ibegin + istep + if (cache_key < X(idx)) then + iend = idx + istep = ishft(idx-ibegin,-1) + cycle + else if (cache_key > X(idx)) then + ibegin = idx + istep = ishft(iend-idx,-1) + cycle + else + value = Y(idx) + return + endif + else if (cache_key > X(idx)) then + ibegin = idx + istep = ishft(iend-idx,-1) + idx = idx + istep + if (cache_key < X(idx)) then + iend = idx + istep = ishft(idx-ibegin,-1) + cycle + else if (cache_key > X(idx)) then + ibegin = idx + istep = ishft(iend-idx,-1) + cycle + else + value = Y(idx) + return + endif + else + value = Y(idx) + return + endif + enddo + idx = ibegin + value = Y(idx) + if (min(iend_in,sze) > ibegin+16) then + iend = ibegin+16 + !DIR$ VECTOR ALIGNED + do while (cache_key > X(idx)) + idx = idx+1 + value = Y(idx) + end do + else + !DIR$ VECTOR ALIGNED + do while (cache_key > X(idx)) + idx = idx+1 + value = Y(idx) + if (idx /= iend) then + cycle + else + exit + endif + end do + endif + if (cache_key /= X(idx)) then + idx = 1-idx + value = 0.d0 + endif + return + + else + + if (cache_key < X(ibegin)) then + idx = -ibegin + value = 0.d0 + return + endif + if (cache_key > X(iend)) then + idx = -iend + value = 0.d0 + return + endif + if (cache_key == X(ibegin)) then + idx = ibegin + value = Y(idx) + return + endif + if (cache_key == X(iend)) then + idx = iend + value = Y(idx) + return + endif + endif + +end + subroutine get_cache_map_n_elements_max(map,n_elements_max) use map_module From 985090e88590a6c3befedd04cceed1db70e1a390 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 25 Nov 2015 10:18:34 +0100 Subject: [PATCH 05/31] Static array sizes in integrals --- src/Integrals_Bielec/ao_bi_integrals.irp.f | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/src/Integrals_Bielec/ao_bi_integrals.irp.f b/src/Integrals_Bielec/ao_bi_integrals.irp.f index 19f4d2a7..34833f3b 100644 --- a/src/Integrals_Bielec/ao_bi_integrals.irp.f +++ b/src/Integrals_Bielec/ao_bi_integrals.irp.f @@ -727,6 +727,7 @@ subroutine integrale_new(I_f,a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z,p,q implicit none + include 'Utils/constants.include.F' double precision :: p,q integer :: a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z integer :: i, n_iter, n_pt, j @@ -741,8 +742,9 @@ subroutine integrale_new(I_f,a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z,p,q p01_1 = 0.5d0/q p10_2 = 0.5d0 * q /(p * q + p * p) p01_2 = 0.5d0 * p /(q * q + q * p) - double precision :: B10(n_pt), B01(n_pt), B00(n_pt) - double precision :: t1(n_pt), t2(n_pt) + double precision :: B00(n_pt_max_integrals) + double precision :: B10(n_pt_max_integrals), B01(n_pt_max_integrals) + double precision :: t1(n_pt_max_integrals), t2(n_pt_max_integrals) !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: t1, t2, B10, B01, B00 ix = a_x+b_x jx = c_x+d_x @@ -797,10 +799,11 @@ recursive subroutine I_x1_new(a,c,B_10,B_01,B_00,res,n_pt) ! recursive function involved in the bielectronic integral END_DOC implicit none + include 'Utils/constants.include.F' integer, intent(in) :: a,c,n_pt - double precision, intent(in) :: B_10(n_pt),B_01(n_pt),B_00(n_pt) - double precision, intent(out) :: res(n_pt) - double precision :: res2(n_pt) + double precision, intent(in) :: B_10(n_pt_max_integrals),B_01(n_pt_max_integrals),B_00(n_pt_max_integrals) + double precision, intent(out) :: res(n_pt_max_integrals) + double precision :: res2(n_pt_max_integrals) integer :: i if(c<0)then @@ -832,9 +835,10 @@ recursive subroutine I_x2_new(c,B_10,B_01,B_00,res,n_pt) BEGIN_DOC ! recursive function involved in the bielectronic integral END_DOC + include 'Utils/constants.include.F' integer, intent(in) :: c, n_pt - double precision, intent(in) :: B_10(n_pt),B_01(n_pt),B_00(n_pt) - double precision, intent(out) :: res(n_pt) + double precision, intent(in) :: B_10(n_pt_max_integrals),B_01(n_pt_max_integrals),B_00(n_pt_max_integrals) + double precision, intent(out) :: res(n_pt_max_integrals) integer :: i if(c==1)then From f0303feb7f6b258272de40d42e91493975c2ff5d Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 25 Nov 2015 10:36:10 +0100 Subject: [PATCH 06/31] Added EZFIO.cfg of MP2 module --- plugins/MP2/EZFIO.cfg | 5 +++++ 1 file changed, 5 insertions(+) create mode 100644 plugins/MP2/EZFIO.cfg diff --git a/plugins/MP2/EZFIO.cfg b/plugins/MP2/EZFIO.cfg new file mode 100644 index 00000000..8577d8f1 --- /dev/null +++ b/plugins/MP2/EZFIO.cfg @@ -0,0 +1,5 @@ +[energy] +type: double precision +doc: MP2 energy +interface: ezfio + From aaa780358e6d436efcd3c41bc252acf2002ce665 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 25 Nov 2015 10:46:53 +0100 Subject: [PATCH 07/31] Typo --- ocaml/qp_run.ml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ocaml/qp_run.ml b/ocaml/qp_run.ml index 2abea107..eb1445d8 100644 --- a/ocaml/qp_run.ml +++ b/ocaml/qp_run.ml @@ -42,7 +42,7 @@ let run exe ezfio_file = let spec = let open Command.Spec in empty - +> anon ("exectuable" %: string) + +> anon ("executable" %: string) +> anon ("ezfio_file" %: string) ;; From 8ca2815964ca24195c836dd7503827517bde8952 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 25 Nov 2015 11:49:42 +0100 Subject: [PATCH 08/31] Added CAS+S --- configure | 30 ++++++--- install/scripts/install_f77zmq.sh | 19 ++++++ install/scripts/install_zeromq.sh | 23 +++++++ plugins/CAS_SD/H_apply.irp.f | 18 +++++- plugins/CAS_SD/cas_s.irp.f | 95 +++++++++++++++++++++++++++++ plugins/CAS_SD/cas_s_selected.irp.f | 89 +++++++++++++++++++++++++++ 6 files changed, 264 insertions(+), 10 deletions(-) create mode 100755 install/scripts/install_f77zmq.sh create mode 100755 install/scripts/install_zeromq.sh create mode 100644 plugins/CAS_SD/cas_s.irp.f create mode 100644 plugins/CAS_SD/cas_s_selected.irp.f diff --git a/configure b/configure index 29c99d24..8376940d 100755 --- a/configure +++ b/configure @@ -63,6 +63,8 @@ d_dependency = { "emsl": ["python"], "gcc": [], "g++": [], + "zeromq" : [ "g++" ], + "f77zmq" : [ "zeromq", "python" ], "python": [], "ninja": ["g++", "python"], "make": [], @@ -121,8 +123,7 @@ ninja = Info( default_path=join(QP_ROOT_BIN, "ninja")) emsl = Info( - url='{head}/LCPQ/EMSL_Basis_Set_Exchange_Local/{tail}'.format(** - path_github), + url='{head}/LCPQ/EMSL_Basis_Set_Exchange_Local/{tail}'.format(**path_github), description=' EMSL basis set library', default_path=join(QP_ROOT_INSTALL, "emsl")) @@ -131,6 +132,16 @@ ezfio = Info( description=' EZFIO', default_path=join(QP_ROOT_INSTALL, "EZFIO")) +zeromq = Info( + url='http://download.zeromq.org/zeromq-4.1.3.tar.gz', + description=' ZeroMQ', + default_path=join(QP_ROOT_INSTALL, "zeromq")) + +f77zmq = Info( + url='{head}/zeromq/f77zmq/{tail}'.format(**path_github), + description=' F77-ZeroMQ', + default_path=join(QP_ROOT_INSTALL, "f77zmq")) + p_graphviz = Info( url='https://github.com/xflr6/graphviz/archive/master.tar.gz', description=' Python library for graphviz', @@ -139,7 +150,8 @@ p_graphviz = Info( d_info = dict() for m in ["ocaml", "m4", "curl", "zlib", "path", "irpf90", "docopt", - "resultsFile", "ninja", "emsl", "ezfio", "p_graphviz"]: + "resultsFile", "ninja", "emsl", "ezfio", "p_graphviz", + "zeromq", "f77zmq" ]: exec ("d_info['{0}']={0}".format(m)) @@ -190,8 +202,7 @@ def check_output(*popenargs, **kwargs): def checking(d_dependency): """ - For each key in d_dependency check if it - is avalabie or not + For each key in d_dependency check if it is avalabie """ def check_python(): @@ -261,7 +272,7 @@ def checking(d_dependency): l_installed = dict() l_needed = [] - # Check all the other + # Check all the others length = max(map(len, d_dependency)) for i in d_dependency.keys(): @@ -276,7 +287,7 @@ def checking(d_dependency): l_needed.append(i) print "" - # Expend the need_stuff for all the genealogy + # Expand the needed stuff for all the genealogy l_install_descendant = get_list_descendant(d_dependency, l_installed, l_needed) @@ -329,7 +340,7 @@ _|_ | | _> |_ (_| | | (_| |_ | (_) | | d_print = { "install_ninja": "Install ninja...", "build": "Creating build.ninja...", - "install": "Installing the dependencies with Ninja..." + "install": "Installing the dependencies using Ninja..." } length = max(map(len, d_print.values())) @@ -373,7 +384,7 @@ _|_ | | _> |_ (_| | | (_| |_ | (_) | | descr = d_info[need].description default_path = d_info[need].default_path - # Build to dowload + # Build to download l_build += ["build {0}: download".format(archive_path), " url = {0}".format(url), " descr = {0}".format(descr), ""] @@ -527,3 +538,4 @@ if __name__ == '__main__': create_ninja_and_rc(l_installed) recommendation() + diff --git a/install/scripts/install_f77zmq.sh b/install/scripts/install_f77zmq.sh new file mode 100755 index 00000000..afbd03ca --- /dev/null +++ b/install/scripts/install_f77zmq.sh @@ -0,0 +1,19 @@ +#!/bin/bash -x + +TARGET=f77zmq + +function _install() +{ + set -e + set -u + export C_INCLUDE_PATH="${C_INCLUDE_PATH}":"${QP_ROOT}"/lib + export ZMQ_H=$PWD/zmq.h + cd "${BUILD}" + make -j 8 + mv libf77zmq.a "${QP_ROOT}"/lib + mv libf77zmq.so "${QP_ROOT}"/lib + cd - + return 0 +} + +source scripts/build.sh diff --git a/install/scripts/install_zeromq.sh b/install/scripts/install_zeromq.sh new file mode 100755 index 00000000..0e36d46c --- /dev/null +++ b/install/scripts/install_zeromq.sh @@ -0,0 +1,23 @@ +#!/bin/bash -x + +TARGET=zeromq + +function _install() +{ + export C_INCLUDE_PATH="${C_INCLUDE_PATH}":./ + set -e + set -u + cd "${BUILD}" + ./configure --without-libsodium || exit 1 + make -j 8 || exit 1 + rm -f -- ../lib/libzmq.a ../lib/libzmq.so ../lib/libzmq.so.5 + cp .libs/libzmq.a "${QP_ROOT}"/lib + cp .libs/libzmq.so "${QP_ROOT}"/lib/libzmq.so.5 + cp include/{zmq.h,zmq_utils.h} "${QP_ROOT}"/lib + cd "${QP_ROOT}"/lib + ln -s libzmq.so.5 libzmq.so + cd - + return 0 +} + +source scripts/build.sh diff --git a/plugins/CAS_SD/H_apply.irp.f b/plugins/CAS_SD/H_apply.irp.f index e2f939fe..35c45fb6 100644 --- a/plugins/CAS_SD/H_apply.irp.f +++ b/plugins/CAS_SD/H_apply.irp.f @@ -5,7 +5,6 @@ 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() @@ -19,5 +18,22 @@ 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/CAS_SD/cas_s.irp.f b/plugins/CAS_SD/cas_s.irp.f new file mode 100644 index 00000000..e0c4a663 --- /dev/null +++ b/plugins/CAS_SD/cas_s.irp.f @@ -0,0 +1,95 @@ +program full_ci + implicit none + integer :: i,k + integer :: N_det_old + + + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) + integer :: N_st, degree + N_st = N_states + allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st)) + character*(64) :: perturbation + PROVIDE N_det_cas + + N_det_old = 0 + pt2 = 1.d0 + diag_algorithm = "Lapack" + 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 + + do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max) + N_det_old = N_det + call H_apply_CAS_S(pt2, norm_pert, H_pert_diag, N_st) + + 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 *, 'N_det = ', N_det + print *, 'N_states = ', N_states + print *, 'PT2 = ', pt2 + print *, 'E = ', CI_energy + print *, 'E+PT2 = ', CI_energy+pt2 + print *, '-----' + call ezfio_set_cas_sd_energy(CI_energy(1)) + if (abort_all) then + exit + endif + if (N_det == N_det_old) then + exit + endif + enddo + call diagonalize_CI + + if(do_pt2_end)then + print*,'Last iteration only to compute the PT2' + threshold_selectors = 1.d0 + threshold_generators = 0.999d0 + call H_apply_CAS_S_PT2(pt2, norm_pert, H_pert_diag, N_st) + + print *, 'Final step' + print *, 'N_det = ', N_det + print *, 'N_states = ', N_states + print *, 'PT2 = ', pt2 + print *, 'E = ', CI_energy + print *, 'E+PT2 = ', CI_energy+pt2 + print *, '-----' + call ezfio_set_cas_sd_energy_pt2(CI_energy(1)+pt2(1)) + endif + + + integer :: exc_max, degree_min + exc_max = 0 + print *, 'CAS determinants : ', N_det_cas + do i=1,min(N_det_cas,10) + do k=i,N_det_cas + call get_excitation_degree(psi_cas(1,1,k),psi_cas(1,1,i),degree,N_int) + exc_max = max(exc_max,degree) + enddo + call debug_det(psi_cas(1,1,i),N_int) + print *, '' + enddo + print *, 'Max excitation degree in the CAS :', exc_max +end diff --git a/plugins/CAS_SD/cas_s_selected.irp.f b/plugins/CAS_SD/cas_s_selected.irp.f new file mode 100644 index 00000000..7a72a243 --- /dev/null +++ b/plugins/CAS_SD/cas_s_selected.irp.f @@ -0,0 +1,89 @@ +program full_ci + implicit none + integer :: i,k + + + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) + integer :: N_st, degree + N_st = N_states + allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st)) + character*(64) :: perturbation + PROVIDE N_det_cas + + pt2 = 1.d0 + diag_algorithm = "Lapack" + 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 + + do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max) + call H_apply_CAS_S_selected(pt2, norm_pert, H_pert_diag, N_st) + + 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 *, 'N_det = ', N_det + print *, 'N_states = ', N_states + print *, 'PT2 = ', pt2 + print *, 'E = ', CI_energy + print *, 'E+PT2 = ', CI_energy+pt2 + print *, '-----' + call ezfio_set_cas_sd_energy(CI_energy(1)) + if (abort_all) then + exit + endif + enddo + call diagonalize_CI + + if(do_pt2_end)then + print*,'Last iteration only to compute the PT2' + threshold_selectors = 1.d0 + threshold_generators = 0.999d0 + call H_apply_CAS_S_PT2(pt2, norm_pert, H_pert_diag, N_st) + + print *, 'Final step' + print *, 'N_det = ', N_det + print *, 'N_states = ', N_states + print *, 'PT2 = ', pt2 + print *, 'E = ', CI_energy + print *, 'E+PT2 = ', CI_energy+pt2 + print *, '-----' + call ezfio_set_cas_sd_energy_pt2(CI_energy(1)+pt2(1)) + endif + + + integer :: exc_max, degree_min + exc_max = 0 + print *, 'CAS determinants : ', N_det_cas + do i=1,min(N_det_cas,10) + do k=i,N_det_cas + call get_excitation_degree(psi_cas(1,1,k),psi_cas(1,1,i),degree,N_int) + exc_max = max(exc_max,degree) + enddo + call debug_det(psi_cas(1,1,i),N_int) + print *, '' + enddo + print *, 'Max excitation degree in the CAS :', exc_max +end From c6dd986ad2f1cd7f7b9da847a61811a518ea10dc Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 25 Nov 2015 13:57:40 +0100 Subject: [PATCH 09/31] ZeroMQ installation OK --- configure | 19 ++++++++++--------- install/scripts/install_f77zmq.sh | 10 +++++----- install/scripts/install_zeromq.sh | 5 +++-- scripts/compilation/qp_create_ninja.py | 3 ++- 4 files changed, 20 insertions(+), 17 deletions(-) diff --git a/configure b/configure index 8376940d..c0925a78 100755 --- a/configure +++ b/configure @@ -46,6 +46,7 @@ if len(sys.argv) != 3: QP_ROOT = os.getcwd() QP_ROOT_BIN = join(QP_ROOT, "bin") +QP_ROOT_LIB = join(QP_ROOT, "lib") QP_ROOT_INSTALL = join(QP_ROOT, "install") os.environ["PATH"] = os.environ["PATH"] + ":" + QP_ROOT_BIN @@ -95,12 +96,12 @@ curl = Info( zlib = Info( url='http://zlib.net/zlib-1.2.8.tar.gz', description=' zlib', - default_path=join(QP_ROOT_INSTALL, "zlib")) + default_path=join(QP_ROOT_LIB, "libz.a")) -path = Info( +patch = Info( url='ftp://ftp.gnu.org/gnu/patch/patch-2.7.5.tar.gz', - description=' path', - default_path=join(QP_ROOT, "lib", "libz.a")) + description=' patch', + default_path=join(QP_ROOT_BIN, "patch")) irpf90 = Info( url='{head}/LCPQ/irpf90/{tail}'.format(**path_github), @@ -119,7 +120,7 @@ resultsFile = Info( ninja = Info( url='{head}/martine/ninja/{tail}'.format(**path_github), - description=' nina', + description=' ninja', default_path=join(QP_ROOT_BIN, "ninja")) emsl = Info( @@ -135,12 +136,12 @@ ezfio = Info( zeromq = Info( url='http://download.zeromq.org/zeromq-4.1.3.tar.gz', description=' ZeroMQ', - default_path=join(QP_ROOT_INSTALL, "zeromq")) + default_path=join(QP_ROOT_LIB, "libzmq.a")) f77zmq = Info( - url='{head}/zeromq/f77zmq/{tail}'.format(**path_github), + url='{head}/zeromq/f77_zmq/{tail}'.format(**path_github), description=' F77-ZeroMQ', - default_path=join(QP_ROOT_INSTALL, "f77zmq")) + default_path=join(QP_ROOT_LIB, "libf77zmq.a")) p_graphviz = Info( url='https://github.com/xflr6/graphviz/archive/master.tar.gz', @@ -149,7 +150,7 @@ p_graphviz = Info( d_info = dict() -for m in ["ocaml", "m4", "curl", "zlib", "path", "irpf90", "docopt", +for m in ["ocaml", "m4", "curl", "zlib", "patch", "irpf90", "docopt", "resultsFile", "ninja", "emsl", "ezfio", "p_graphviz", "zeromq", "f77zmq" ]: exec ("d_info['{0}']={0}".format(m)) diff --git a/install/scripts/install_f77zmq.sh b/install/scripts/install_f77zmq.sh index afbd03ca..e8d592e1 100755 --- a/install/scripts/install_f77zmq.sh +++ b/install/scripts/install_f77zmq.sh @@ -4,14 +4,14 @@ TARGET=f77zmq function _install() { + export C_INCLUDE_PATH="${C_INCLUDE_PATH}":"${QP_ROOT}"/lib set -e set -u - export C_INCLUDE_PATH="${C_INCLUDE_PATH}":"${QP_ROOT}"/lib - export ZMQ_H=$PWD/zmq.h + export ZMQ_H="${QP_ROOT}"/lib/zmq.h cd "${BUILD}" - make -j 8 - mv libf77zmq.a "${QP_ROOT}"/lib - mv libf77zmq.so "${QP_ROOT}"/lib + make -j 8 || exit 1 + mv libf77zmq.a "${QP_ROOT}"/lib || exit 1 + mv libf77zmq.so "${QP_ROOT}"/lib || exit 1 cd - return 0 } diff --git a/install/scripts/install_zeromq.sh b/install/scripts/install_zeromq.sh index 0e36d46c..5ea280a7 100755 --- a/install/scripts/install_zeromq.sh +++ b/install/scripts/install_zeromq.sh @@ -7,16 +7,17 @@ function _install() export C_INCLUDE_PATH="${C_INCLUDE_PATH}":./ set -e set -u + ORIG=$(pwd) cd "${BUILD}" ./configure --without-libsodium || exit 1 make -j 8 || exit 1 - rm -f -- ../lib/libzmq.a ../lib/libzmq.so ../lib/libzmq.so.5 + rm -f -- "${QP_ROOT}"/lib/libzmq.a "${QP_ROOT}"/lib/libzmq.so "${QP_ROOT}"/lib/libzmq.so.5 cp .libs/libzmq.a "${QP_ROOT}"/lib cp .libs/libzmq.so "${QP_ROOT}"/lib/libzmq.so.5 cp include/{zmq.h,zmq_utils.h} "${QP_ROOT}"/lib cd "${QP_ROOT}"/lib ln -s libzmq.so.5 libzmq.so - cd - + cd ${ORIG} return 0 } diff --git a/scripts/compilation/qp_create_ninja.py b/scripts/compilation/qp_create_ninja.py index b0fe2f8c..7171c2df 100755 --- a/scripts/compilation/qp_create_ninja.py +++ b/scripts/compilation/qp_create_ninja.py @@ -37,6 +37,7 @@ 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") ROOT_BUILD_NINJA = join(QP_ROOT, "config", "build.ninja") header = r"""# @@ -95,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") - l_string.append("LIB = {0} {1} {2}".format(LIB, lib_lapack, EZFIO_LIB)) + l_string.append("LIB = {0} {1} {2} {3}".format(LIB, lib_lapack, EZFIO_LIB, ZMQ_LIB)) l_string.append("") From 6cad30269f44890ba3506009098475b7caa5f409 Mon Sep 17 00:00:00 2001 From: Michel Caffarel Date: Wed, 25 Nov 2015 18:04:05 +0100 Subject: [PATCH 10/31] Guess based on the overlap matrix of AOs to remove linear dependencies --- src/MOGuess/guess_overlap.irp.f | 14 ++++++++++++++ src/MOGuess/truncate_mos.irp.f | 10 ++++++++++ src/MO_Basis/utils.irp.f | 23 +++++++++++++++++++++++ 3 files changed, 47 insertions(+) create mode 100644 src/MOGuess/guess_overlap.irp.f create mode 100644 src/MOGuess/truncate_mos.irp.f diff --git a/src/MOGuess/guess_overlap.irp.f b/src/MOGuess/guess_overlap.irp.f new file mode 100644 index 00000000..768724c6 --- /dev/null +++ b/src/MOGuess/guess_overlap.irp.f @@ -0,0 +1,14 @@ +program guess_mimi + BEGIN_DOC +! Produce `H_core` MO orbital + END_DOC + implicit none + character*(64) :: label + mo_coef = ao_ortho_lowdin_coef + TOUCH mo_coef + label = "Guess" + call mo_as_eigvectors_of_mo_matrix(-ao_overlap, & + size(ao_overlap,1), & + size(ao_overlap,2),label) + call save_mos +end diff --git a/src/MOGuess/truncate_mos.irp.f b/src/MOGuess/truncate_mos.irp.f new file mode 100644 index 00000000..29756055 --- /dev/null +++ b/src/MOGuess/truncate_mos.irp.f @@ -0,0 +1,10 @@ +program prog_truncate_mo + BEGIN_DOC +! Truncate MO set + END_DOC + implicit none + integer :: n + write(*,*) 'Number of MOs to keep' + read (*,*) n + call save_mos_truncated(n) +end diff --git a/src/MO_Basis/utils.irp.f b/src/MO_Basis/utils.irp.f index 7cc94c6d..6f96ab93 100644 --- a/src/MO_Basis/utils.irp.f +++ b/src/MO_Basis/utils.irp.f @@ -21,6 +21,29 @@ subroutine save_mos end +subroutine save_mos_truncated(n) + implicit none + double precision, allocatable :: buffer(:,:) + integer :: i,j,n + + call system('$QP_ROOT/scripts/save_current_mos.sh '//trim(ezfio_filename)) + + call ezfio_set_mo_basis_mo_tot_num(n) + call ezfio_set_mo_basis_mo_label(mo_label) + call ezfio_set_mo_basis_ao_md5(ao_md5) + allocate ( buffer(ao_num,n) ) + buffer = 0.d0 + do j = 1, n + do i = 1, ao_num + buffer(i,j) = mo_coef(i,j) + enddo + enddo + call ezfio_set_mo_basis_mo_coef(buffer) + call ezfio_set_mo_basis_mo_occ(mo_occ) + deallocate (buffer) + +end + subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label) implicit none integer,intent(in) :: n,m From 109e3946eb046fb2a7431d6df022c9f8bb15b697 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 25 Nov 2015 18:39:49 +0100 Subject: [PATCH 11/31] Added print in configure for debugging --- install/scripts/build.sh | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/install/scripts/build.sh b/install/scripts/build.sh index 6b7fc80a..835cdc11 100755 --- a/install/scripts/build.sh +++ b/install/scripts/build.sh @@ -5,6 +5,11 @@ BUILD=_build/${TARGET} rm -rf -- ${BUILD} mkdir ${BUILD} || exit 1 tar -zxf Downloads/${TARGET}.tar.gz --strip-components=1 --directory=${BUILD} || exit 1 -_install || exit 1 +_install +if [[ $? -ne 0 ]] +then + cat _build/${TARGET}.log + exit 1 +fi rm -rf -- ${BUILD} _build/${TARGET}.log -exit 0 \ No newline at end of file +exit 0 From 6b02b622ba4c502c9a9b7c4ed9e61eb6cc5ef2c8 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 25 Nov 2015 18:55:06 +0100 Subject: [PATCH 12/31] Added print in configure for debugging --- configure | 11 ++++++++++- install/scripts/build.sh | 7 +------ 2 files changed, 11 insertions(+), 7 deletions(-) diff --git a/configure b/configure index c0925a78..84291ad8 100755 --- a/configure +++ b/configure @@ -417,7 +417,16 @@ _|_ | | _> |_ (_| | | (_| |_ | (_) | | path_ninja = find_path("ninja", l_installed) subprocess.check_call("cd install ;{0}".format(path_ninja), shell=True) except: - raise + prefix = os.path.join('install', '_build') + for filename in os.listdir(prefix): + if filename.endswith(".log"): + with open( os.path.join(prefix,filename) ,'r') as f: + print "\n\n" + print "=-=-=-=-=-=- %s =-=-=-=-=-=-" %(filename) + print f.read() + print "=-=-=-=-=-=-=-=-=-=-=-=-=-=-\n\n" + print "Error in installation of dependencies" + sys.exit(1) else: print r""" _________ diff --git a/install/scripts/build.sh b/install/scripts/build.sh index 835cdc11..79a71065 100755 --- a/install/scripts/build.sh +++ b/install/scripts/build.sh @@ -5,11 +5,6 @@ BUILD=_build/${TARGET} rm -rf -- ${BUILD} mkdir ${BUILD} || exit 1 tar -zxf Downloads/${TARGET}.tar.gz --strip-components=1 --directory=${BUILD} || exit 1 -_install -if [[ $? -ne 0 ]] -then - cat _build/${TARGET}.log - exit 1 -fi +_install || exit 1 rm -rf -- ${BUILD} _build/${TARGET}.log exit 0 From d3530303ed72f96fb334319180392ee0a880cb67 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 25 Nov 2015 18:59:49 +0100 Subject: [PATCH 13/31] QP_ROOT in install scripts --- install/scripts/install_f77zmq.sh | 3 +++ install/scripts/install_zeromq.sh | 3 +++ 2 files changed, 6 insertions(+) diff --git a/install/scripts/install_f77zmq.sh b/install/scripts/install_f77zmq.sh index e8d592e1..ff2417fb 100755 --- a/install/scripts/install_f77zmq.sh +++ b/install/scripts/install_f77zmq.sh @@ -4,6 +4,9 @@ TARGET=f77zmq function _install() { + cd .. + QP_ROOT=$PWD + cd - export C_INCLUDE_PATH="${C_INCLUDE_PATH}":"${QP_ROOT}"/lib set -e set -u diff --git a/install/scripts/install_zeromq.sh b/install/scripts/install_zeromq.sh index 5ea280a7..9508f457 100755 --- a/install/scripts/install_zeromq.sh +++ b/install/scripts/install_zeromq.sh @@ -4,6 +4,9 @@ TARGET=zeromq function _install() { + cd .. + QP_ROOT=$PWD + cd - export C_INCLUDE_PATH="${C_INCLUDE_PATH}":./ set -e set -u From 89a376512a766bab9ba7b4a2428dbb7cd03df35e Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 25 Nov 2015 23:46:09 +0100 Subject: [PATCH 14/31] Added list of basis sets --- ocaml/qp_create_ezfio_from_xyz.ml | 28 ++++++++++++++++++++++++---- 1 file changed, 24 insertions(+), 4 deletions(-) diff --git a/ocaml/qp_create_ezfio_from_xyz.ml b/ocaml/qp_create_ezfio_from_xyz.ml index b2eaaae0..ee7e7d40 100644 --- a/ocaml/qp_create_ezfio_from_xyz.ml +++ b/ocaml/qp_create_ezfio_from_xyz.ml @@ -17,7 +17,7 @@ let spec = ~doc:"int Spin multiplicity (2S+1) of the molecule. Default is 1." +> flag "p" no_arg ~doc:" Using pseudopotentials" - +> anon ("xyz_file" %: string) + +> anon ("xyz_file" %: file ) let dummy_centers ~threshold ~molecule ~nuclei = @@ -68,8 +68,21 @@ let dummy_centers ~threshold ~molecule ~nuclei = ) - - +let list_basis () = + let basis_list = + Qpackage.root ^ "/install/emsl/EMSL_api.py list_basis" + |> Unix.open_process_in + |> In_channel.input_lines + |> List.map ~f:(fun x -> + match String.split x ~on:'\'' with + | [] -> "" + | a :: [] + | _ :: a :: _ -> String.strip a + ) + in + List.sort basis_list ~cmp:String.ascending + |> String.concat ~sep:"\t" + let run ?o b c d m p xyz_file = @@ -345,6 +358,13 @@ let command = Command.basic ~summary: "Quantum Package command" ~readme:(fun () -> " + +=== Available basis sets === + +" ^ (list_basis ()) ^ " + +============================ + Creates an EZFIO directory from a standard xyz file. The basis set is defined as a single string if all the atoms are taken from the same basis set, otherwise specific elements can be defined as follows: @@ -354,7 +374,7 @@ otherwise specific elements can be defined as follows: If a file with the same name as the basis set exists, this file will be read. Otherwise, the basis set is obtained from the database. - ") +" ) spec (fun o b c d m p xyz_file () -> run ?o b c d m p xyz_file ) From 47ac033293b0cda5281dd94dd584c141a1d7f695 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 26 Nov 2015 14:54:40 +0100 Subject: [PATCH 15/31] Removed dead code in pseudopotentials --- src/Integrals_Monoelec/pseudopot.f90 | 140 --------------------------- 1 file changed, 140 deletions(-) diff --git a/src/Integrals_Monoelec/pseudopot.f90 b/src/Integrals_Monoelec/pseudopot.f90 index c262105f..32402c74 100644 --- a/src/Integrals_Monoelec/pseudopot.f90 +++ b/src/Integrals_Monoelec/pseudopot.f90 @@ -1585,148 +1585,8 @@ end bessel_mod_exp=x**n*accu end -! double precision function bessel_mod(x,n) -! IMPLICIT DOUBLE PRECISION (A-H,O-Z) -! parameter(NBESSMAX=100) -! dimension SI(0:NBESSMAX),DI(0:NBESSMAX) -! if(n.lt.0.or.n.gt.NBESSMAX)stop 'pb with argument of bessel_mod' -! CALL SPHI(N,X,NBESSMAX,SI,DI) -! bessel_mod=si(n) -! end - - SUBROUTINE SPHI(N,X,NMAX,SI,DI) -! -! ======================================================== -! Purpose: Compute modified spherical Bessel functions -! of the first kind, in(x) and in'(x) -! Input : x --- Argument of in(x) -! n --- Order of in(x) ( n = 0,1,2,... ) -! Output: SI(n) --- in(x) -! DI(n) --- in'(x) -! NM --- Highest order computed -! Routines called: -! MSTA1 and MSTA2 for computing the starting -! point for backward recurrence -! ======================================================== -! - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DIMENSION SI(0:NMAX),DI(0:NMAX) - NM=N - IF (DABS(X).LT.1.0D-100) THEN - DO 10 K=0,N - SI(K)=0.0D0 -10 DI(K)=0.0D0 - SI(0)=1.0D0 - DI(1)=0.333333333333333D0 - RETURN - ENDIF - SI(0)=DSINH(X)/X - SI(1)=-(DSINH(X)/X-DCOSH(X))/X - SI0=SI(0) - IF (N.GE.2) THEN - M=MSTA1(X,200) - - write(34,*)'m=',m - - IF (M.LT.N) THEN - NM=M - ELSE - M=MSTA2(X,N,15) - write(34,*)'m=',m - ENDIF - F0=0.0D0 - F1=1.0D0-100 - DO 15 K=M,0,-1 - F=(2.0D0*K+3.0D0)*F1/X+F0 - IF (K.LE.NM) SI(K)=F - F0=F1 -15 F1=F - CS=SI0/F - write(34,*)'cs=',cs - DO 20 K=0,NM -20 SI(K)=CS*SI(K) - ENDIF - DI(0)=SI(1) - DO 25 K=1,NM -25 DI(K)=SI(K-1)-(K+1.0D0)/X*SI(K) - RETURN - END - INTEGER FUNCTION MSTA1(X,MP) -! -! =================================================== -! Purpose: Determine the starting point for backward -! recurrence such that the magnitude of -! Jn(x) at that point is about 10^(-MP) -! Input : x --- Argument of Jn(x) -! MP --- Value of magnitude -! Output: MSTA1 --- Starting point -! =================================================== -! - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - A0=DABS(X) - N0=INT(1.1*A0)+1 - F0=ENVJ(N0,A0)-MP - N1=N0+5 - F1=ENVJ(N1,A0)-MP - DO 10 IT=1,20 - NN=N1-(N1-N0)/(1.0D0-F0/F1) - F=ENVJ(NN,A0)-MP - IF(ABS(NN-N1).LT.1) GO TO 20 - N0=N1 - F0=F1 - N1=NN - 10 F1=F - 20 MSTA1=NN - RETURN - END - - - INTEGER FUNCTION MSTA2(X,N,MP) -! -! =================================================== -! Purpose: Determine the starting point for backward -! recurrence such that all Jn(x) has MP -! significant digits -! Input : x --- Argument of Jn(x) -! n --- Order of Jn(x) -! MP --- Significant digit -! Output: MSTA2 --- Starting point -! =================================================== -! - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - A0=DABS(X) - HMP=0.5D0*MP - EJN=ENVJ(N,A0) - IF (EJN.LE.HMP) THEN - OBJ=MP - N0=INT(1.1*A0) - ELSE - OBJ=HMP+EJN - N0=N - ENDIF - F0=ENVJ(N0,A0)-OBJ - N1=N0+5 - F1=ENVJ(N1,A0)-OBJ - DO 10 IT=1,20 - NN=N1-(N1-N0)/(1.0D0-F0/F1) - F=ENVJ(NN,A0)-OBJ - IF (iABS(NN-N1).LT.1) GO TO 20 - N0=N1 - F0=F1 - N1=NN -10 F1=F -20 MSTA2=NN+10 - RETURN - END - - double precision FUNCTION ENVJ(N,X) - DOUBLE PRECISION X - integer N - ENVJ=0.5D0*DLOG10(6.28D0*N)-N*DLOG10(1.36D0*X/N) - RETURN - END !c Computation of real spherical harmonics Ylm(theta,phi) !c From 61bcd8c10328bfc0ab94909dc2be5b3dfb3f805c Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 26 Nov 2015 15:00:02 +0100 Subject: [PATCH 16/31] Eigenvalues are now printed with the correct sign --- plugins/Hartree_Fock/damping_SCF.irp.f | 2 +- plugins/Hartree_Fock/huckel.irp.f | 2 +- src/Determinants/density_matrix.irp.f | 4 +-- src/MOGuess/H_CORE_guess.irp.f | 2 +- src/MOGuess/guess_overlap.irp.f | 4 +-- src/MOGuess/h_core_guess_routine.irp.f | 2 +- src/MO_Basis/utils.irp.f | 37 +++++++++++++++++++------- 7 files changed, 36 insertions(+), 17 deletions(-) diff --git a/plugins/Hartree_Fock/damping_SCF.irp.f b/plugins/Hartree_Fock/damping_SCF.irp.f index c7c005c9..d77c91c5 100644 --- a/plugins/Hartree_Fock/damping_SCF.irp.f +++ b/plugins/Hartree_Fock/damping_SCF.irp.f @@ -119,7 +119,7 @@ subroutine damping_SCF write(output_hartree_fock,'(A4,X,A16, X, A16, X, A16, X, A4 )'), '====','================','================','================', '====' write(output_hartree_fock,*) - call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label) + call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1) call write_double(output_hartree_fock, E_min, 'Hartree-Fock energy') call ezfio_set_hartree_fock_energy(E_min) diff --git a/plugins/Hartree_Fock/huckel.irp.f b/plugins/Hartree_Fock/huckel.irp.f index 6139ac46..8f61f0cf 100644 --- a/plugins/Hartree_Fock/huckel.irp.f +++ b/plugins/Hartree_Fock/huckel.irp.f @@ -13,7 +13,7 @@ subroutine huckel_guess label = "Guess" call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral, & size(mo_mono_elec_integral,1), & - size(mo_mono_elec_integral,2),label) + size(mo_mono_elec_integral,2),label,1) TOUCH mo_coef c = 0.5d0 * 1.75d0 diff --git a/src/Determinants/density_matrix.irp.f b/src/Determinants/density_matrix.irp.f index 5f087642..9aeb658e 100644 --- a/src/Determinants/density_matrix.irp.f +++ b/src/Determinants/density_matrix.irp.f @@ -185,9 +185,9 @@ subroutine set_natural_mos allocate(tmp(size(one_body_dm_mo,1),size(one_body_dm_mo,2))) ! Negation to have the occupied MOs first after the diagonalization - tmp = -one_body_dm_mo + tmp = one_body_dm_mo label = "Natural" - call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label) + call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,-1) deallocate(tmp) end diff --git a/src/MOGuess/H_CORE_guess.irp.f b/src/MOGuess/H_CORE_guess.irp.f index 1893c08b..b65fe07d 100644 --- a/src/MOGuess/H_CORE_guess.irp.f +++ b/src/MOGuess/H_CORE_guess.irp.f @@ -10,7 +10,7 @@ program H_CORE_guess label = "Guess" call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral, & size(mo_mono_elec_integral,1), & - size(mo_mono_elec_integral,2),label) + size(mo_mono_elec_integral,2),label,1) print *, 'save mos' call save_mos diff --git a/src/MOGuess/guess_overlap.irp.f b/src/MOGuess/guess_overlap.irp.f index 768724c6..cf222507 100644 --- a/src/MOGuess/guess_overlap.irp.f +++ b/src/MOGuess/guess_overlap.irp.f @@ -7,8 +7,8 @@ program guess_mimi mo_coef = ao_ortho_lowdin_coef TOUCH mo_coef label = "Guess" - call mo_as_eigvectors_of_mo_matrix(-ao_overlap, & + call mo_as_eigvectors_of_mo_matrix(ao_overlap, & size(ao_overlap,1), & - size(ao_overlap,2),label) + size(ao_overlap,2),label,-1) call save_mos end diff --git a/src/MOGuess/h_core_guess_routine.irp.f b/src/MOGuess/h_core_guess_routine.irp.f index 566592ba..605c7c8a 100644 --- a/src/MOGuess/h_core_guess_routine.irp.f +++ b/src/MOGuess/h_core_guess_routine.irp.f @@ -9,7 +9,7 @@ subroutine hcore_guess label = "Guess" call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral, & size(mo_mono_elec_integral,1), & - size(mo_mono_elec_integral,2),label) + size(mo_mono_elec_integral,2),label,1) print *, 'save mos' call save_mos SOFT_TOUCH mo_coef mo_label diff --git a/src/MO_Basis/utils.irp.f b/src/MO_Basis/utils.irp.f index 6f96ab93..3f4a260a 100644 --- a/src/MO_Basis/utils.irp.f +++ b/src/MO_Basis/utils.irp.f @@ -44,13 +44,14 @@ subroutine save_mos_truncated(n) end -subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label) +subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label,sign) implicit none - integer,intent(in) :: n,m + integer,intent(in) :: n,m, sign character*(64), intent(in) :: label double precision, intent(in) :: matrix(n,m) - double precision, allocatable :: mo_coef_new(:,:), R(:,:),eigvalues(:) + integer :: i,j + double precision, allocatable :: mo_coef_new(:,:), R(:,:),eigvalues(:), A(:,:) !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: mo_coef_new, R call write_time(output_mo_basis) @@ -58,30 +59,48 @@ subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label) print *, irp_here, ': Error : m/= mo_tot_num' stop 1 endif - allocate(R(n,m),mo_coef_new(ao_num_align,m),eigvalues(m)) + allocate(A(n,m),R(n,m),mo_coef_new(ao_num_align,m),eigvalues(m)) + if (sign == -1) then + do j=1,m + do i=1,n + A(i,j) = -matrix(i,j) + enddo + enddo + else + do j=1,m + do i=1,n + A(i,j) = matrix(i,j) + enddo + enddo + endif mo_coef_new = mo_coef - call lapack_diag(eigvalues,R,matrix,size(matrix,1),size(matrix,2)) - integer :: i + call lapack_diag(eigvalues,R,A,n,m) write (output_mo_basis,'(A)'), 'MOs are now **'//trim(label)//'**' write (output_mo_basis,'(A)'), '' write (output_mo_basis,'(A)'), 'Eigenvalues' write (output_mo_basis,'(A)'), '-----------' write (output_mo_basis,'(A)'), '' write (output_mo_basis,'(A)'), '======== ================' - do i = 1, m - write (output_mo_basis,'(I8,X,F16.10)'), i,eigvalues(i) + if (sign == -1) then + do i=1,m + eigvalues(i) = -eigvalues(i) + enddo + endif + do i=1,m + write (output_mo_basis,'(I8,X,F16.10)'), i,eigvalues(i) enddo write (output_mo_basis,'(A)'), '======== ================' write (output_mo_basis,'(A)'), '' call dgemm('N','N',ao_num,m,m,1.d0,mo_coef_new,size(mo_coef_new,1),R,size(R,1),0.d0,mo_coef,size(mo_coef,1)) - deallocate(mo_coef_new,R,eigvalues) + deallocate(A,mo_coef_new,R,eigvalues) 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) implicit none integer,intent(in) :: n,m From 8b664eb8f382a1d32814a4b312b5037880f4219f Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 27 Nov 2015 09:02:20 +0100 Subject: [PATCH 17/31] qp_install_module creates a default main and installs --- scripts/module/qp_install_module.py | 55 ++++++++++++++++++++++++++--- 1 file changed, 50 insertions(+), 5 deletions(-) diff --git a/scripts/module/qp_install_module.py b/scripts/module/qp_install_module.py index 2dddfcb7..b656b95a 100755 --- a/scripts/module/qp_install_module.py +++ b/scripts/module/qp_install_module.py @@ -59,10 +59,48 @@ def save_new_module(path, l_child): f.write(D_KEY["needed_module"]) f.write(D_KEY["documentation"]) + with open(os.path.join(path, "%s.main.irp.f"%(module_name) ), "w") as f: + f.write("program {0}".format(module_name) ) + f.write(""" implicit none + BEGIN_DOC +! TODO + END_DOC + print *, ' _/ ' + print *, ' -:\_?, _Jm####La ' + print *, 'J"(:" > _]#AZ#Z#UUZ##, ' + print *, '_,::./ %(|i%12XmX1*1XL _?, ' + print *, ' \..\ _\(vmWQwodY+ia%lnL _",/ ( ' + print *, ' .:< ]J=mQD?WXn|,)nr" ' + print *, ' 4XZ#Xov1v}=)vnXAX1nnv;1n" ' + print *, ' ]XX#ZXoovvvivnnnlvvo2*i7 ' + print *, ' "23Z#1S2oo2XXSnnnoSo2>v" ' + print *, ' miX#L -~`""!!1}oSoe|i7 ' + print *, ' 4cn#m, v221=|v[ ' + print *, ' ]hI3Zma,;..__wXSe=+vo ' + print *, ' ]Zov*XSUXXZXZXSe||vo2 ' + print *, ' ]Z#>=|< ' + print *, ' -ziiiii||||||+||==+> ' + print *, ' -%|+++||=|=+|=|==/ ' + print *, ' -a>====+|====-:- ' + print *, ' "~,- -- /- ' + print *, ' -. )> ' + print *, ' .~ +- ' + print *, ' . .... : . ' + print *, ' -------~ ' + print *, '' +end +""") -if __name__ == '__main__': - arguments = docopt(__doc__) - +def main(arguments): if arguments["list"]: if arguments["--installed"]: @@ -109,12 +147,14 @@ if __name__ == '__main__': save_new_module(path, l_child_reduce) print " [ OK ]" - print "Your module is created in the `plugins` directory." - print "You need to create some `.irp.f` to be able to install it." # print "` {0} install {1} `".format(os.path.basename(__file__), name) print "" + arguments["create"]=False + arguments["install"]=True + main(arguments) elif arguments["download"]: + print "Not yet implemented" pass # d_local = get_dict_child([QP_SRC]) # d_remote = get_dict_child(arguments[""]) @@ -207,3 +247,8 @@ if __name__ == '__main__': except OSError: print "%s is a core module which can't be removed" % module + +if __name__ == '__main__': + arguments = docopt(__doc__) + main(arguments) + From 627d405e54b5283a186df9ebb3fdaeea074aea49 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 27 Nov 2015 09:03:03 +0100 Subject: [PATCH 18/31] Default level shift is 0.5 au --- plugins/Hartree_Fock/EZFIO.cfg | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/plugins/Hartree_Fock/EZFIO.cfg b/plugins/Hartree_Fock/EZFIO.cfg index 1b8486ee..d8207cc4 100644 --- a/plugins/Hartree_Fock/EZFIO.cfg +++ b/plugins/Hartree_Fock/EZFIO.cfg @@ -14,7 +14,7 @@ default: 200 type: Positive_float doc: Energy shift on the virtual MOs to improve SCF convergence interface: ezfio,provider,ocaml -default: 0.0 +default: 0.5 [mo_guess_type] type: MO_guess From 281c6aae3343779daad7967d9ab4fde2329efade Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 27 Nov 2015 09:20:28 +0100 Subject: [PATCH 19/31] Added tput after vim --- ocaml/.gitignore | 1 - ocaml/qp_edit.ml | 308 +++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 308 insertions(+), 1 deletion(-) create mode 100644 ocaml/qp_edit.ml diff --git a/ocaml/.gitignore b/ocaml/.gitignore index 5618a6c0..651376e2 100644 --- a/ocaml/.gitignore +++ b/ocaml/.gitignore @@ -13,7 +13,6 @@ qp_basis_clean.native qp_create_ezfio_from_xyz qp_create_ezfio_from_xyz.native qp_edit -qp_edit.ml qp_edit.native qp_print qp_print.native diff --git a/ocaml/qp_edit.ml b/ocaml/qp_edit.ml new file mode 100644 index 00000000..0ad55ef0 --- /dev/null +++ b/ocaml/qp_edit.ml @@ -0,0 +1,308 @@ +open Qputils;; +open Qptypes;; +open Core.Std;; + +(** Interactive editing of the input. + +WARNING +This file is autogenerad by +`${QP_ROOT}/script/ezfio_interface/ei_handler.py` +*) + + +(** Keywords used to define input sections *) +type keyword = +| Ao_basis +| Determinants_by_hand +| Electrons +| Mo_basis +| Nuclei +| Determinants +| Hartree_fock +| Integrals_bielec +| Perturbation +| Properties +| Pseudo +;; + + +let keyword_to_string = function +| Ao_basis -> "AO basis" +| Determinants_by_hand -> "Determinants_by_hand" +| Electrons -> "Electrons" +| Mo_basis -> "MO basis" +| Nuclei -> "Molecule" +| Determinants -> "Determinants" +| Hartree_fock -> "Hartree_fock" +| Integrals_bielec -> "Integrals_bielec" +| Perturbation -> "Perturbation" +| Properties -> "Properties" +| Pseudo -> "Pseudo" +;; + + + +(** Create the header of the temporary file *) +let file_header filename = + Printf.sprintf " +================================================================== + Quantum Package +================================================================== + +Editing file `%s` + +" filename +;; + + +(** Creates the header of a section *) +let make_header kw = + let s = keyword_to_string kw in + let l = String.length s in + "\n\n"^s^"\n"^(String.init l ~f:(fun _ -> '='))^"\n\n" +;; + + +(** Returns the rst string of section [s] *) +let get s = + let header = (make_header s) in + let f (read,to_rst) = + match read () with + | Some text -> header ^ (Rst_string.to_string (to_rst text)) + | None -> "" + in + let rst = + try + begin + let open Input in + match s with + | Mo_basis -> + f Mo_basis.(read, to_rst) + | Electrons -> + f Electrons.(read, to_rst) + | Nuclei -> + f Nuclei.(read, to_rst) + | Ao_basis -> + f Ao_basis.(read, to_rst) + | Determinants_by_hand -> + f Determinants_by_hand.(read, to_rst) + | Determinants -> + f Determinants.(read, to_rst) + | Hartree_fock -> + f Hartree_fock.(read, to_rst) + | Integrals_bielec -> + f Integrals_bielec.(read, to_rst) + | Perturbation -> + f Perturbation.(read, to_rst) + | Properties -> + f Properties.(read, to_rst) + | Pseudo -> + f Pseudo.(read, to_rst) + end + with + | Sys_error msg -> (Printf.eprintf "Info: %s\n%!" msg ; "") + in + rst +;; + + +(** Applies the changes from the string [str] corresponding to section [s] *) +let set str s = + let header = (make_header s) in + match String.substr_index ~pos:0 ~pattern:header str with + | None -> () + | Some idx -> + begin + let index_begin = idx + (String.length header) in + let index_end = + match ( String.substr_index ~pos:(index_begin+(String.length header)+1) + ~pattern:"==" str) with + | Some i -> i + | None -> String.length str + in + let l = index_end - index_begin in + let str = String.sub ~pos:index_begin ~len:l str + |> Rst_string.of_string + in + let write (of_rst,w) s = + try + match of_rst str with + | Some data -> w data + | None -> () + with + | _ -> (Printf.eprintf "Info: Read error in %s\n%!" + (keyword_to_string s); ignore (of_rst str) ) + in + let open Input in + match s with + | Determinants -> write Determinants.(of_rst, write) s + | Hartree_fock -> write Hartree_fock.(of_rst, write) s + | Integrals_bielec -> write Integrals_bielec.(of_rst, write) s + | Perturbation -> write Perturbation.(of_rst, write) s + | Properties -> write Properties.(of_rst, write) s + | Pseudo -> write Pseudo.(of_rst, write) s + | Electrons -> write Electrons.(of_rst, write) s + | Determinants_by_hand -> write Determinants_by_hand.(of_rst, write) s + | Nuclei -> write Nuclei.(of_rst, write) s + | Ao_basis -> () (* TODO *) + | Mo_basis -> () (* TODO *) + end +;; + + +(** Creates the temporary file for interactive editing *) +let create_temp_file ezfio_filename fields = + let temp_filename = Filename.temp_file "qp_edit_" ".rst" in + begin + Out_channel.with_file temp_filename ~f:(fun out_channel -> + (file_header ezfio_filename) :: (List.map ~f:get fields) + |> String.concat ~sep:"\n" + |> Out_channel.output_string out_channel + ) + end + ; temp_filename +;; + + + +let run check_only ezfio_filename = + + (* Open EZFIO *) + if (not (Sys.file_exists_exn ezfio_filename)) then + failwith (ezfio_filename^" does not exists"); + + Ezfio.set_file ezfio_filename; + + (* + let output = (file_header ezfio_filename) :: ( + List.map ~f:get [ + Ao_basis ; + Mo_basis ; + ]) + in + String.concat output + |> print_string + *) + + let tasks = [ + Nuclei ; + Ao_basis; + Electrons ; + Determinants ; + Hartree_fock ; + Integrals_bielec ; + Perturbation ; + Properties ; + Pseudo ; + Mo_basis; + Determinants_by_hand ; + ] + in + + (* Create the temp file *) + let temp_filename = create_temp_file ezfio_filename tasks in + + (* Open the temp file with external editor *) + let editor = + match Sys.getenv "EDITOR" with + | Some editor -> editor + | None -> "vi" + in + + match check_only with + | true -> () + | false -> + Printf.sprintf "%s %s ; tput init &> /dev/null" editor temp_filename + |> Sys.command_exn + ; + + (* Re-read the temp file *) + let temp_string = + In_channel.with_file temp_filename ~f:(fun in_channel -> + In_channel.input_all in_channel) + in + List.iter ~f:(fun x -> set temp_string x) tasks; + + (* Remove temp_file *) + Sys.remove temp_filename; +;; + + +(** Create a backup file in case of an exception *) +let create_backup ezfio_filename = + Printf.sprintf " + rm -f %s/backup.tgz ; + tar -zcf .backup.tgz %s && mv .backup.tgz %s/backup.tgz + " + ezfio_filename ezfio_filename ezfio_filename + |> Sys.command_exn +;; + + +(** Restore the backup file when an exception occuprs *) +let restore_backup ezfio_filename = + Printf.sprintf "tar -zxf %s/backup.tgz" + ezfio_filename + |> Sys.command_exn +;; + + +let spec = + let open Command.Spec in + empty + +> flag "-c" no_arg + ~doc:"Checks the input data" +(* + +> flag "o" (optional string) + ~doc:"Prints output data" +*) + +> anon ("ezfio_file" %: string) +;; + +let command = + Command.basic + ~summary: "Quantum Package command" + ~readme:(fun () -> + " +Edit input data + ") + spec +(* (fun i o ezfio_file () -> *) + (*fun ezfio_file () -> + try + run ezfio_file + with + | _ msg -> print_string ("\n\nError\n\n"^msg^"\n\n") + *) + (fun c ezfio_file () -> + try + run c ezfio_file ; + (* create_backup ezfio_file; *) + with + | Failure exc + | Invalid_argument exc as e -> + begin + Printf.eprintf "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-\n\n"; + Printf.eprintf "%s\n\n" exc; + Printf.eprintf "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-\n\n"; + (* restore_backup ezfio_file; *) + raise e + end + | Assert_failure (file, line, ch) as e -> + begin + Printf.eprintf "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-\n\n"; + Printf.eprintf "Assert error in file $QP_ROOT/ocaml/%s, line %d, character %d\n\n" file line ch; + Printf.eprintf "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-\n\n"; + (* restore_backup ezfio_file; *) + raise e + end + ) +;; + +let () = + Command.run command; + exit 0 +;; + + + From 4855837019e800e987151829588492b39bc4354c Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 27 Nov 2015 09:37:16 +0100 Subject: [PATCH 20/31] tput init -> tput sgr0 --- ocaml/qp_edit.ml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ocaml/qp_edit.ml b/ocaml/qp_edit.ml index 0ad55ef0..f6a2ac9c 100644 --- a/ocaml/qp_edit.ml +++ b/ocaml/qp_edit.ml @@ -212,7 +212,7 @@ let run check_only ezfio_filename = match check_only with | true -> () | false -> - Printf.sprintf "%s %s ; tput init &> /dev/null" editor temp_filename + Printf.sprintf "%s %s ; tput sgr0 2> /dev/null" editor temp_filename |> Sys.command_exn ; From 7f4634e49b983897dd4332c5f2f6abca455d046d Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 27 Nov 2015 10:15:46 +0100 Subject: [PATCH 21/31] Isolated SVD routine --- ocaml/.gitignore | 1 + src/AO_Basis/ao_overlap.irp.f | 2 ++ src/Utils/LinearAlgebra.irp.f | 66 +++++++++++++++++++++++++---------- 3 files changed, 50 insertions(+), 19 deletions(-) diff --git a/ocaml/.gitignore b/ocaml/.gitignore index 651376e2..5618a6c0 100644 --- a/ocaml/.gitignore +++ b/ocaml/.gitignore @@ -13,6 +13,7 @@ qp_basis_clean.native qp_create_ezfio_from_xyz qp_create_ezfio_from_xyz.native qp_edit +qp_edit.ml qp_edit.native qp_print qp_print.native diff --git a/src/AO_Basis/ao_overlap.irp.f b/src/AO_Basis/ao_overlap.irp.f index 737f03f7..9b3fa7e9 100644 --- a/src/AO_Basis/ao_overlap.irp.f +++ b/src/AO_Basis/ao_overlap.irp.f @@ -59,6 +59,8 @@ enddo enddo !$OMP END PARALLEL DO + + ! Check for linear dependencies in the basis set END_PROVIDER diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f index a7462f94..c5e66006 100644 --- a/src/Utils/LinearAlgebra.irp.f +++ b/src/Utils/LinearAlgebra.irp.f @@ -1,3 +1,48 @@ +subroutine svd(A,LDA,U,LDU,D,Vt,LDVt,m,n) + implicit none + BEGIN_DOC + ! Compute A = U.D.Vt + ! + ! LDx : leftmost dimension of x + ! + ! Dimsneion of A is m x n + ! + END_DOC + + integer, intent(in) :: LDA, LDU, LDVt, m, n + double precision, intent(in) :: A(LDA,n) + double precision, intent(out) :: U(LDU,n) + double precision,intent(out) :: Vt(LDVt,n) + double precision,intent(out) :: D(n) + double precision,allocatable :: work(:) + integer :: info, lwork, i, j, k + + double precision,allocatable :: A_tmp(:,:) + allocate (A_tmp(LDA,n)) + A_tmp = A + + ! Find optimal size for temp arrays + allocate(work(1)) + lwork = -1 + call dgesvd('A','A', n, n, A_tmp, LDA, & + D, U, LDU, Vt, LDVt, work, lwork, info) + lwork = work(1) + deallocate(work) + + allocate(work(lwork)) + call dgesvd('A','A', n, n, A_tmp, LDA, & + D, U, LDU, Vt, LDVt, work, lwork, info) + deallocate(work,A_tmp) + + if (info /= 0) then + print *, info, ': SVD failed' + stop + endif + +end + + + subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m) implicit none BEGIN_DOC @@ -29,25 +74,8 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m) !DEC$ ATTRIBUTES ALIGN : 64 :: U, Vt, D, work integer :: info, lwork, i, j, k - double precision,allocatable :: overlap_tmp(:,:) - allocate (overlap_tmp(lda,n)) - overlap_tmp = overlap - - allocate(work(1)) - lwork = -1 - call dgesvd('A','A', n, n, overlap_tmp, lda, & - D, U, ldc, Vt, lda, work, lwork, info) - lwork = work(1) - deallocate(work) - allocate(work(lwork)) - call dgesvd('A','A', n, n, overlap_tmp, lda, & - D, U, ldc, Vt, lda, work, lwork, info) - deallocate(work,overlap_tmp) - if (info /= 0) then - print *, info, ': SVD failed' - stop - endif - + call svd(overlap,lda,U,ldc,D,Vt,lda,m,n) + !$OMP PARALLEL DEFAULT(NONE) & !$OMP SHARED(S_half,U,D,Vt,n,C,m) & !$OMP PRIVATE(i,j,k) From 88c47b1f73370cdad670b16807c32d6ff8d7c6d8 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 27 Nov 2015 10:30:50 +0100 Subject: [PATCH 22/31] Trap linear dependencies in basis set --- src/AO_Basis/ao_overlap.irp.f | 1 - src/MOGuess/guess_overlap.irp.f | 1 + src/Utils/LinearAlgebra.irp.f | 4 ++-- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/AO_Basis/ao_overlap.irp.f b/src/AO_Basis/ao_overlap.irp.f index 9b3fa7e9..4487ff77 100644 --- a/src/AO_Basis/ao_overlap.irp.f +++ b/src/AO_Basis/ao_overlap.irp.f @@ -60,7 +60,6 @@ enddo !$OMP END PARALLEL DO - ! Check for linear dependencies in the basis set END_PROVIDER diff --git a/src/MOGuess/guess_overlap.irp.f b/src/MOGuess/guess_overlap.irp.f index cf222507..c2f090e5 100644 --- a/src/MOGuess/guess_overlap.irp.f +++ b/src/MOGuess/guess_overlap.irp.f @@ -4,6 +4,7 @@ program guess_mimi END_DOC implicit none character*(64) :: label + mo_coef = ao_ortho_lowdin_coef TOUCH mo_coef label = "Guess" diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f index c5e66006..dab9e921 100644 --- a/src/Utils/LinearAlgebra.irp.f +++ b/src/Utils/LinearAlgebra.irp.f @@ -82,8 +82,8 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m) !$OMP DO do i=1,n - if ( D(i) < 1.d-6 ) then - D(i) = 0.d0 + if ( D(i) < 1.d-12 ) then + stop 'Linear dependence in basis set' else D(i) = 1.d0/dsqrt(D(i)) endif From ad4dc469cc156a870f258eb0d9dffb5490e31f9d Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 27 Nov 2015 11:42:46 +0100 Subject: [PATCH 23/31] SCF works with non-square matrices --- plugins/Hartree_Fock/Fock_matrix.irp.f | 47 ++++++++++++++++--- .../Hartree_Fock/HF_density_matrix_ao.irp.f | 4 +- 2 files changed, 42 insertions(+), 9 deletions(-) diff --git a/plugins/Hartree_Fock/Fock_matrix.irp.f b/plugins/Hartree_Fock/Fock_matrix.irp.f index 82540de3..12ee276b 100644 --- a/plugins/Hartree_Fock/Fock_matrix.irp.f +++ b/plugins/Hartree_Fock/Fock_matrix.irp.f @@ -344,30 +344,47 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_ao, (ao_num_align, ao_num) ] enddo enddo else - double precision, allocatable :: T(:,:), M(:,:) + double precision, allocatable :: T(:,:), M(:,:) + integer :: ierr ! F_ao = S C F_mo C^t S - allocate (T(ao_num_align,ao_num),M(ao_num_align,ao_num)) - call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & + allocate (T(ao_num_align,ao_num),M(ao_num_align,ao_num),stat=ierr) + if (ierr /=0 ) then + print *, irp_here, ' : allocation failed' + endif + +! ao_overlap (ao_num,ao_num) . mo_coef (ao_num,mo_tot_num) +! -> M(ao_num,mo_tot_num) + call dgemm('N','N', ao_num,mo_tot_num,ao_num, 1.d0, & ao_overlap, size(ao_overlap,1), & mo_coef, size(mo_coef,1), & 0.d0, & M, size(M,1)) + +! M(ao_num,mo_tot_num) . Fock_matrix_mo (mo_tot_num,mo_tot_num) +! -> T(ao_num,mo_tot_num) call dgemm('N','N', ao_num,mo_tot_num,mo_tot_num, 1.d0, & M, size(M,1), & Fock_matrix_mo, size(Fock_matrix_mo,1), & 0.d0, & T, size(T,1)) - call dgemm('N','T', mo_tot_num,ao_num,mo_tot_num, 1.d0, & + +! T(ao_num,mo_tot_num) . mo_coef^T (mo_tot_num,ao_num) +! -> M(ao_num,ao_num) + call dgemm('N','T', ao_num,ao_num,mo_tot_num, 1.d0, & T, size(T,1), & mo_coef, size(mo_coef,1), & 0.d0, & M, size(M,1)) + +! M(ao_num,ao_num) . ao_overlap (ao_num,ao_num) +! -> Fock_matrix_ao(ao_num,ao_num) call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & M, size(M,1), & ao_overlap, size(ao_overlap,1), & 0.d0, & Fock_matrix_ao, size(Fock_matrix_ao,1)) + deallocate(T) endif END_PROVIDER @@ -380,23 +397,39 @@ subroutine Fock_mo_to_ao(FMO,LDFMO,FAO,LDFAO) double precision, intent(out) :: FAO(LDFAO,*) double precision, allocatable :: T(:,:), M(:,:) + integer :: ierr ! F_ao = S C F_mo C^t S - allocate (T(ao_num_align,ao_num),M(ao_num_align,ao_num)) - call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & + allocate (T(ao_num_align,ao_num),M(ao_num_align,ao_num),stat=ierr) + if (ierr /=0 ) then + print *, irp_here, ' : allocation failed' + endif + +! ao_overlap (ao_num,ao_num) . mo_coef (ao_num,mo_tot_num) +! -> M(ao_num,mo_tot_num) + call dgemm('N','N', ao_num,mo_tot_num,ao_num, 1.d0, & ao_overlap, size(ao_overlap,1), & mo_coef, size(mo_coef,1), & 0.d0, & M, size(M,1)) + +! M(ao_num,mo_tot_num) . FMO (mo_tot_num,mo_tot_num) +! -> T(ao_num,mo_tot_num) call dgemm('N','N', ao_num,mo_tot_num,mo_tot_num, 1.d0, & M, size(M,1), & FMO, size(FMO,1), & 0.d0, & T, size(T,1)) - call dgemm('N','T', mo_tot_num,ao_num,mo_tot_num, 1.d0, & + +! T(ao_num,mo_tot_num) . mo_coef^T (mo_tot_num,ao_num) +! -> M(ao_num,ao_num) + call dgemm('N','T', ao_num,ao_num,mo_tot_num, 1.d0, & T, size(T,1), & mo_coef, size(mo_coef,1), & 0.d0, & M, size(M,1)) + +! M(ao_num,ao_num) . ao_overlap (ao_num,ao_num) +! -> Fock_matrix_ao(ao_num,ao_num) call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & M, size(M,1), & ao_overlap, size(ao_overlap,1), & diff --git a/plugins/Hartree_Fock/HF_density_matrix_ao.irp.f b/plugins/Hartree_Fock/HF_density_matrix_ao.irp.f index e8585f59..597b88c7 100644 --- a/plugins/Hartree_Fock/HF_density_matrix_ao.irp.f +++ b/plugins/Hartree_Fock/HF_density_matrix_ao.irp.f @@ -4,7 +4,7 @@ BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_alpha, (ao_num_align,ao_ ! S^-1 x Alpha density matrix in the AO basis x S^-1 END_DOC - call dgemm('N','T',ao_num,ao_num,elec_alpha_num,1.d0, & + call dgemm('N','T',ao_num,mo_tot_num,elec_alpha_num,1.d0, & mo_coef, size(mo_coef,1), & mo_coef, size(mo_coef,1), 0.d0, & HF_density_matrix_ao_alpha, size(HF_density_matrix_ao_alpha,1)) @@ -17,7 +17,7 @@ BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_beta, (ao_num_align,ao_ ! S^-1 Beta density matrix in the AO basis x S^-1 END_DOC - call dgemm('N','T',ao_num,ao_num,elec_beta_num,1.d0, & + call dgemm('N','T',ao_num,mo_tot_num,elec_beta_num,1.d0, & mo_coef, size(mo_coef,1), & mo_coef, size(mo_coef,1), 0.d0, & HF_density_matrix_ao_beta, size(HF_density_matrix_ao_beta,1)) From e1e2a1d3e446d3a9898675eea3de734bb6ae7328 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 27 Nov 2015 15:20:15 +0100 Subject: [PATCH 24/31] Corrected bug for large PT2 --- plugins/Perturbation/perturbation.template.f | 24 +++++++++++++++----- 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/plugins/Perturbation/perturbation.template.f b/plugins/Perturbation/perturbation.template.f index b41c7685..33bd10dd 100644 --- a/plugins/Perturbation/perturbation.template.f +++ b/plugins/Perturbation/perturbation.template.f @@ -20,14 +20,18 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c integer, external :: connected_to_ref logical, external :: is_in_wavefunction - integer(bit_kind) :: minilist(Nint,2,N_det_selectors) - integer :: idx_minilist(N_det_selectors), N_minilist + integer(bit_kind), allocatable :: minilist(:,:,:) + integer, allocatable :: idx_minilist(:) + integer :: N_minilist - integer(bit_kind) :: minilist_gen(Nint,2,N_det_generators) + integer(bit_kind), allocatable :: minilist_gen(:,:,:) integer :: N_minilist_gen logical :: fullMatch logical, external :: is_connected_to + allocate( minilist(Nint,2,N_det_selectors), & + minilist_gen(Nint,2,N_det_generators), & + idx_minilist(N_det_selectors) ) ASSERT (Nint > 0) @@ -40,6 +44,7 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c call create_minilist_find_previous(key_mask, psi_det_generators, miniList_gen, i_generator-1, N_minilist_gen, fullMatch, Nint) if(fullMatch) then + deallocate( minilist, minilist_gen, idx_minilist ) return end if @@ -66,6 +71,7 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c enddo enddo + deallocate( minilist, minilist_gen, idx_minilist ) end @@ -89,14 +95,18 @@ subroutine perturb_buffer_by_mono_$PERT(i_generator,buffer,buffer_size,e_2_pert_ integer, external :: connected_to_ref_by_mono logical, external :: is_in_wavefunction - integer(bit_kind) :: minilist(Nint,2,N_det_selectors) - integer :: idx_minilist(N_det_selectors), N_minilist + integer(bit_kind), allocatable :: minilist(:,:,:) + integer, allocatable :: idx_minilist(:) + integer :: N_minilist - integer(bit_kind) :: minilist_gen(Nint,2,N_det_generators) + integer(bit_kind), allocatable :: minilist_gen(:,:,:) integer :: N_minilist_gen logical :: fullMatch logical, external :: is_connected_to + allocate( minilist(Nint,2,N_det_selectors), & + minilist_gen(Nint,2,N_det_generators), & + idx_minilist(N_det_selectors) ) ASSERT (Nint > 0) @@ -109,6 +119,7 @@ subroutine perturb_buffer_by_mono_$PERT(i_generator,buffer,buffer_size,e_2_pert_ call create_minilist_find_previous(key_mask, psi_det_generators, miniList_gen, i_generator-1, N_minilist_gen, fullMatch, Nint) if(fullMatch) then + deallocate( minilist, minilist_gen, idx_minilist ) return end if @@ -137,6 +148,7 @@ subroutine perturb_buffer_by_mono_$PERT(i_generator,buffer,buffer_size,e_2_pert_ enddo enddo + deallocate( minilist, minilist_gen, idx_minilist ) end From dd78dd43248f5b1195e571f39d2fb7a50feffa5b Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 27 Nov 2015 19:14:36 +0100 Subject: [PATCH 25/31] Hartree-Fock with less MOs than AOs OK --- .../Hartree_Fock/HF_density_matrix_ao.irp.f | 4 +- plugins/Hartree_Fock/diagonalize_fock.irp.f | 129 ++++++++++++------ 2 files changed, 89 insertions(+), 44 deletions(-) diff --git a/plugins/Hartree_Fock/HF_density_matrix_ao.irp.f b/plugins/Hartree_Fock/HF_density_matrix_ao.irp.f index 597b88c7..e8585f59 100644 --- a/plugins/Hartree_Fock/HF_density_matrix_ao.irp.f +++ b/plugins/Hartree_Fock/HF_density_matrix_ao.irp.f @@ -4,7 +4,7 @@ BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_alpha, (ao_num_align,ao_ ! S^-1 x Alpha density matrix in the AO basis x S^-1 END_DOC - call dgemm('N','T',ao_num,mo_tot_num,elec_alpha_num,1.d0, & + call dgemm('N','T',ao_num,ao_num,elec_alpha_num,1.d0, & mo_coef, size(mo_coef,1), & mo_coef, size(mo_coef,1), 0.d0, & HF_density_matrix_ao_alpha, size(HF_density_matrix_ao_alpha,1)) @@ -17,7 +17,7 @@ BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_beta, (ao_num_align,ao_ ! S^-1 Beta density matrix in the AO basis x S^-1 END_DOC - call dgemm('N','T',ao_num,mo_tot_num,elec_beta_num,1.d0, & + call dgemm('N','T',ao_num,ao_num,elec_beta_num,1.d0, & mo_coef, size(mo_coef,1), & mo_coef, size(mo_coef,1), 0.d0, & HF_density_matrix_ao_beta, size(HF_density_matrix_ao_beta,1)) diff --git a/plugins/Hartree_Fock/diagonalize_fock.irp.f b/plugins/Hartree_Fock/diagonalize_fock.irp.f index 7f25851c..850ba0aa 100644 --- a/plugins/Hartree_Fock/diagonalize_fock.irp.f +++ b/plugins/Hartree_Fock/diagonalize_fock.irp.f @@ -10,58 +10,103 @@ integer, allocatable :: iwork(:) double precision, allocatable :: work(:), F(:,:), S(:,:) - allocate(F(ao_num_align,ao_num), S(ao_num_align,ao_num) ) - do j=1,ao_num - do i=1,ao_num - S(i,j) = ao_overlap(i,j) - F(i,j) = Fock_matrix_ao(i,j) - enddo - enddo - n = ao_num - lwork = 1+6*n + 2*n*n - liwork = 3 + 5*n - - allocate(work(lwork), iwork(liwork) ) + if (mo_tot_num == ao_num) then + ! Solve H.C = E.S.C in AO basis set - lwork = -1 - liwork = -1 + allocate(F(ao_num_align,ao_num), S(ao_num_align,ao_num) ) + do j=1,ao_num + do i=1,ao_num + S(i,j) = ao_overlap(i,j) + F(i,j) = Fock_matrix_ao(i,j) + enddo + enddo - call dsygvd(1,'v','u',ao_num,F,size(F,1),S,size(S,1),& - diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info) -! call dsygv(1, 'v', 'u',ao_num,F,size(F,1),S,size(S,1),& -! diagonal_Fock_matrix_mo, work, lwork, info) + n = ao_num + lwork = 1+6*n + 2*n*n + liwork = 3 + 5*n + + allocate(work(lwork), iwork(liwork) ) + lwork = -1 + liwork = -1 + call dsygvd(1,'v','u',ao_num,F,size(F,1),S,size(S,1),& + diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info) - if (info /= 0) then - print *, irp_here//' failed : ', info - stop 1 - endif - lwork = int(work(1)) - liwork = iwork(1) - deallocate(work,iwork) - allocate(work(lwork), iwork(liwork) ) -! deallocate(work) -! allocate(work(lwork)) + if (info /= 0) then + print *, irp_here//' failed : ', info + stop 1 + endif + lwork = int(work(1)) + liwork = iwork(1) + deallocate(work,iwork) + allocate(work(lwork), iwork(liwork) ) - call dsygvd(1,'v','u',ao_num,F,size(F,1),S,size(S,1),& - diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info) + call dsygvd(1,'v','u',ao_num,F,size(F,1),S,size(S,1),& + diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info) -! call dsygv(1, 'v', 'u',ao_num,F,size(F,1),S,size(S,1),& -! diagonal_Fock_matrix_mo, work, lwork, info) + if (info /= 0) then + print *, irp_here//' failed : ', info + stop 1 + endif + do j=1,mo_tot_num + do i=1,ao_num + eigenvectors_Fock_matrix_mo(i,j) = F(i,j) + enddo + enddo - if (info /= 0) then - print *, irp_here//' failed : ', info - stop 1 - endif - do j=1,mo_tot_num - do i=1,ao_num - eigenvectors_Fock_matrix_mo(i,j) = F(i,j) - enddo - enddo + deallocate(work, iwork, F, S) + + else + + ! Solve H.C = E.C in MO basis set + + allocate( F(mo_tot_num_align,mo_tot_num) ) + do j=1,mo_tot_num + do i=1,mo_tot_num + F(i,j) = Fock_matrix_mo(i,j) + enddo + enddo + + n = mo_tot_num + lwork = 1+6*n + 2*n*n + liwork = 3 + 5*n + + allocate(work(lwork), iwork(liwork) ) + + lwork = -1 + liwork = -1 + + call dsyevd( 'V', 'U', mo_tot_num, F, & + size(F,1), diagonal_Fock_matrix_mo, & + work, lwork, iwork, liwork, info) + + if (info /= 0) then + print *, irp_here//' failed : ', info + stop 1 + endif + lwork = int(work(1)) + liwork = iwork(1) + deallocate(work,iwork) + allocate(work(lwork), iwork(liwork) ) + + call dsyevd( 'V', 'U', mo_tot_num, F, & + size(F,1), diagonal_Fock_matrix_mo, & + work, lwork, iwork, liwork, info) + + if (info /= 0) then + print *, irp_here//' failed : ', info + stop 1 + endif + + call dgemm('N','N',ao_num,mo_tot_num,mo_tot_num, 1.d0, & + mo_coef, size(mo_coef,1), F, size(F,1), & + 0.d0, eigenvectors_Fock_matrix_mo, size(eigenvectors_Fock_matrix_mo,1)) + deallocate(work, iwork, F) + + endif - deallocate(work, iwork, F, S) END_PROVIDER BEGIN_PROVIDER [double precision, diagonal_Fock_matrix_mo_sum, (mo_tot_num)] From 1d58caea77f2630b002efa751d7ac08dcb8c6a18 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 30 Nov 2015 14:54:37 +0100 Subject: [PATCH 26/31] 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 27/31] 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 28/31] 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 29/31] 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 30/31] 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 31/31] 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 +