diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index 517220a87..af7420c85 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -178,7 +178,7 @@ subroutine select_singles_and_doubles(i_generator, hole_mask, particle_mask, foc integer(bit_kind), allocatable :: minilist(:, :, :), fullminilist(:, :, :) logical, allocatable :: banned(:,:,:), bannedOrb(:,:) double precision, allocatable :: coef_fullminilist_rev(:,:) - double precision, allocatable :: mat(:,:,:), hij_cache(:,:,:) + double precision, allocatable :: mat(:,:,:) PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique @@ -205,7 +205,7 @@ subroutine select_singles_and_doubles(i_generator, hole_mask, particle_mask, foc ! Removed to avoid introducing determinants already presents in the wf !double precision, parameter :: norm_thr = 1.d-16 - allocate (indices(N_det), hij_cache(mo_num,mo_num,2), & + allocate (indices(N_det), & exc_degree(max(N_det_alpha_unique,N_det_beta_unique))) ! Pre-compute excitation degrees wrt alpha determinants @@ -511,15 +511,10 @@ subroutine select_singles_and_doubles(i_generator, hole_mask, particle_mask, foc maskInd = maskInd + 1 if(mod(maskInd, csubset) == (subset-1)) then - call get_mo_two_e_integrals_ij(h2,h1,mo_num,hij_cache(1,1,1),mo_integrals_map) - if (sp /= 3) then ! AA or BB - call get_mo_two_e_integrals_ij(h1,h2,mo_num,hij_cache(1,1,2),mo_integrals_map) - endif - call spot_isinwf(mask, fullminilist, i_generator, fullinteresting(0), banned, fullMatch, fullinteresting) if(fullMatch) cycle - call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting, hij_cache) + call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting) call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf) end if @@ -535,7 +530,7 @@ subroutine select_singles_and_doubles(i_generator, hole_mask, particle_mask, foc enddo enddo deallocate(preinteresting, prefullinteresting, interesting, fullinteresting) - deallocate(banned, bannedOrb, mat, hij_cache) + deallocate(banned, bannedOrb, mat) end subroutine BEGIN_TEMPLATE @@ -765,6 +760,7 @@ subroutine fill_buffer_$DOUBLE(i_generator, sp, h1, h2, bannedOrb, banned, fock_ enddo do_diag = sum(dabs(coef)) > 0.001d0 .and. N_states > 1 + !do_diag = .False. double precision :: eigvalues(N_states+1) double precision :: work(1+6*(N_states+1)+2*(N_states+1)**2) @@ -918,7 +914,7 @@ subroutine fill_buffer_$DOUBLE(i_generator, sp, h1, h2, bannedOrb, banned, fock_ END_TEMPLATE -subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, interesting, hij_cache) +subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, interesting) use bitmasks implicit none BEGIN_DOC @@ -930,7 +926,6 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere integer, intent(in) :: sp, i_gen, N_sel integer, intent(in) :: interesting(0:N_sel) integer(bit_kind),intent(in) :: mask(N_int, 2), det(N_int, 2, N_sel) - double precision, intent(in) :: hij_cache(mo_num, mo_num, 2) logical, intent(inout) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num, 2) double precision, intent(inout) :: mat(N_states, mo_num, mo_num) @@ -1000,9 +995,9 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere if(nt == 4) then call get_d2(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i))) else if(nt == 3) then - call get_d1(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)), hij_cache) + call get_d1(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i))) else - call get_d0(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)), hij_cache) + call get_d0(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i))) end if else if(nt == 4) then call bitstring_to_list_in_selection(mobMask(1,1), p(1,1), p(0,1), N_int) @@ -1203,7 +1198,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) end -subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs, hij_cache) +subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) use bitmasks implicit none @@ -1214,7 +1209,6 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs, double precision, intent(in) :: coefs(N_states) double precision, intent(inout) :: mat(N_states, mo_num, mo_num) integer, intent(in) :: h(0:2,2), p(0:4,2), sp - double precision, intent(in) :: hij_cache(mo_num, mo_num, 2) double precision, external :: get_phase_bi, mo_two_e_integral logical :: ok @@ -1256,11 +1250,13 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs, p1 = p(1,ma) p2 = p(2,ma) if(.not. bannedOrb(puti, mi)) then + call get_mo_two_e_integrals(hfix,p1,p2,mo_num,hij_cache1(1,1),mo_integrals_map) + call get_mo_two_e_integrals(hfix,p2,p1,mo_num,hij_cache1(1,2),mo_integrals_map) tmp_row = 0d0 do putj=1, hfix-1 if(lbanned(putj, ma)) cycle if(banned(putj, puti,bant)) cycle - hij = hij_cache(hfix,putj,1) - hij_cache(putj,hfix,1) + hij = hij_cache1(putj,1) - hij_cache1(putj,2) if (hij /= 0.d0) then hij = hij * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int) !DIR$ LOOP COUNT AVG(4) @@ -1272,7 +1268,7 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs, do putj=hfix+1, mo_num if(lbanned(putj, ma)) cycle if(banned(putj, puti,bant)) cycle - hij = hij_cache(putj,hfix,1) - hij_cache(hfix,putj,1) + hij = hij_cache1(putj,2) - hij_cache1(putj,1) if (hij /= 0.d0) then hij = hij * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int) !DIR$ LOOP COUNT AVG(4) @@ -1463,7 +1459,7 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs, -subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs, hij_cache) +subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) use bitmasks implicit none @@ -1474,7 +1470,6 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs, double precision, intent(in) :: coefs(N_states) double precision, intent(inout) :: mat(N_states, mo_num, mo_num) integer, intent(in) :: h(0:2,2), p(0:4,2), sp - double precision, intent(in) :: hij_cache(mo_num, mo_num, 2) integer :: i, j, k, s, h1, h2, p1, p2, puti, putj double precision :: hij, phase @@ -1483,6 +1478,9 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs, integer, parameter :: bant=1 + double precision, allocatable :: hij_cache1(:), hij_cache2(:) + allocate (hij_cache1(mo_num),hij_cache2(mo_num)) + PROVIDE mo_integrals_threshold if(sp == 3) then ! AB @@ -1490,6 +1488,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs, h2 = p(1,2) do p1=1, mo_num if(bannedOrb(p1, 1)) cycle + call get_mo_two_e_integrals(p1,h2,h1,mo_num,hij_cache1,mo_integrals_map) do p2=1, mo_num if(bannedOrb(p2,2)) cycle if(banned(p1, p2, bant)) cycle ! rentable? @@ -1498,7 +1497,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs, call i_h_j(gen, det, N_int, hij) else phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int) - hij = hij_cache(p2,p1,1) * phase + hij = hij_cache1(p2) * phase end if if (dabs(hij) < mo_integrals_threshold) cycle !DIR$ LOOP COUNT AVG(4) @@ -1513,6 +1512,8 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs, p2 = p(2,sp) do puti=1, mo_num if (bannedOrb(puti, sp)) cycle + call get_mo_two_e_integrals(puti,p2,p1,mo_num,hij_cache1,mo_integrals_map) + call get_mo_two_e_integrals(puti,p1,p2,mo_num,hij_cache2,mo_integrals_map) do putj=puti+1, mo_num if(bannedOrb(putj, sp)) cycle if(banned(puti, putj, bant)) cycle ! rentable? @@ -1521,7 +1522,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs, call i_h_j(gen, det, N_int, hij) if (dabs(hij) < mo_integrals_threshold) cycle else - hij = hij_cache(putj,puti,1) - hij_cache(putj,puti,2) + hij = hij_cache1(putj) - hij_cache2(putj) if (dabs(hij) < mo_integrals_threshold) cycle hij = hij * get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2, N_int) end if @@ -1533,6 +1534,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs, end do end if + deallocate(hij_cache1,hij_cache2) end diff --git a/src/cipsi_utils/pt2_stoch_routines.irp.f b/src/cipsi_utils/pt2_stoch_routines.irp.f index 100335f69..744c4006f 100644 --- a/src/cipsi_utils/pt2_stoch_routines.irp.f +++ b/src/cipsi_utils/pt2_stoch_routines.irp.f @@ -894,9 +894,8 @@ integer function pt2_find_sample_lr(v, w, l_in, r_in) do t=1, pt2_N_teeth tooth_width = tilde_cW(pt2_n_0(t+1)) - tilde_cW(pt2_n_0(t)) if (tooth_width == 0.d0) then - tooth_width = sum(tilde_w(pt2_n_0(t):pt2_n_0(t+1))) + tooth_width = max(1.d-15,sum(tilde_w(pt2_n_0(t):pt2_n_0(t+1)))) endif - ASSERT(tooth_width > 0.d0) do i=pt2_n_0(t)+1, pt2_n_0(t+1) pt2_w(i) = tilde_w(i) * pt2_W_T / tooth_width end do diff --git a/src/mo_two_e_ints/mo_bi_integrals.irp.f b/src/mo_two_e_ints/mo_bi_integrals.irp.f index 6ee23f4a3..dd5428980 100644 --- a/src/mo_two_e_ints/mo_bi_integrals.irp.f +++ b/src/mo_two_e_ints/mo_bi_integrals.irp.f @@ -70,6 +70,9 @@ subroutine mo_two_e_integrals_index(i,j,k,l,i1) else call add_integrals_to_map(full_ijkl_bitmask_4) endif + + integer*8 :: get_mo_map_size, mo_map_size + mo_map_size = get_mo_map_size() double precision, external :: map_mb print*,'Molecular integrals provided:' print*,' Size of MO map ', map_mb(mo_integrals_map) ,'MB' @@ -79,9 +82,6 @@ subroutine mo_two_e_integrals_index(i,j,k,l,i1) call wall_time(wall_2) call cpu_time(cpu_2) - integer*8 :: get_mo_map_size, mo_map_size - mo_map_size = get_mo_map_size() - print*,' cpu time :',cpu_2 - cpu_1, 's' print*,' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1), ')' diff --git a/src/utils/mmap.f90 b/src/utils/mmap.f90 index 4ac322333..826ea9197 100644 --- a/src/utils/mmap.f90 +++ b/src/utils/mmap.f90 @@ -2,6 +2,8 @@ module mmap_module use iso_c_binding + character*(256) :: mmap_prefix = '/tmp/' + type mmap_type type(c_ptr) :: ptr ! Pointer to the data character*(128) :: filename ! Name of the file @@ -155,7 +157,15 @@ subroutine mmap_create(filename, shape, bytes, read_only, single_node, map) map%filename = filename else call getenv('EZFIO_FILE', map%filename) - map%filename = trim(map%filename) // '/work/tmpfile' + if (trim(map%filename) /= '') then + map%filename = trim(map%filename) // '/work/' + else + call getenv('TMPDIR', map%filename) + if (trim(map%filename) == '') then + map%filename = '/tmp/' + endif + endif + map%filename = trim(map%filename) // '/tmpfile' endif map%length = int(bytes,8)