diff --git a/plugins/local/bi_ort_ints/total_twoe_pot.irp.f b/plugins/local/bi_ort_ints/total_twoe_pot.irp.f index 5e6a24e9..bf5cc36f 100644 --- a/plugins/local/bi_ort_ints/total_twoe_pot.irp.f +++ b/plugins/local/bi_ort_ints/total_twoe_pot.irp.f @@ -40,38 +40,95 @@ BEGIN_PROVIDER [double precision, mo_bi_ortho_tc_two_e_chemist, (mo_num, mo_num, END_DOC implicit none - integer :: i, j, k, l, m, n, p, q + integer :: i, j, k, l, m, n, p, q, s, r + double precision :: t1, t2 double precision, allocatable :: a1(:,:,:,:), a2(:,:,:,:) + double precision, allocatable :: a_jkp(:,:,:), a_kpq(:,:,:), a_pqr(:,:,:) + + print *, ' PROVIDING mo_bi_ortho_tc_two_e_chemist ...' + call wall_time(t1) + call print_memory_usage() PROVIDE mo_r_coef mo_l_coef + PROVIDe ao_two_e_tc_tot - allocate(a2(ao_num,ao_num,ao_num,mo_num)) + if(ao_to_mo_tc_n3) then - call dgemm( 'T', 'N', ao_num*ao_num*ao_num, mo_num, ao_num, 1.d0 & - , ao_two_e_tc_tot(1,1,1,1), ao_num, mo_l_coef(1,1), ao_num & - , 0.d0 , a2(1,1,1,1), ao_num*ao_num*ao_num) + print*, ' memory scale of TC ao -> mo: O(N3) ' - allocate(a1(ao_num,ao_num,mo_num,mo_num)) + allocate(a_jkp(ao_num,ao_num,mo_num)) + allocate(a_kpq(ao_num,mo_num,mo_num)) + allocate(a_pqr(mo_num,mo_num,mo_num)) - call dgemm( 'T', 'N', ao_num*ao_num*mo_num, mo_num, ao_num, 1.d0 & - , a2(1,1,1,1), ao_num, mo_r_coef(1,1), ao_num & - , 0.d0, a1(1,1,1,1), ao_num*ao_num*mo_num) + do s = 1, mo_num + mo_bi_ortho_tc_two_e_chemist(:,:,:,s) = 0.d0 - deallocate(a2) - allocate(a2(ao_num,mo_num,mo_num,mo_num)) + do l = 1, ao_num - call dgemm( 'T', 'N', ao_num*mo_num*mo_num, mo_num, ao_num, 1.d0 & - , a1(1,1,1,1), ao_num, mo_l_coef(1,1), ao_num & - , 0.d0, a2(1,1,1,1), ao_num*mo_num*mo_num) + call dgemm( 'T', 'N', ao_num*ao_num, mo_num, ao_num, 1.d0 & + , ao_two_e_tc_tot(1,1,1,l), ao_num, mo_l_coef(1,1), ao_num & + , 0.d0, a_jkp(1,1,1), ao_num*ao_num) + + call dgemm( 'T', 'N', ao_num*mo_num, mo_num, ao_num, 1.d0 & + , a_jkp(1,1,1), ao_num, mo_r_coef(1,1), ao_num & + , 0.d0, a_kpq(1,1,1), ao_num*mo_num) + + call dgemm( 'T', 'N', mo_num*mo_num, mo_num, ao_num, 1.d0 & + , a_kpq(1,1,1), ao_num, mo_l_coef(1,1), ao_num & + , 0.d0, a_pqr(1,1,1), mo_num*mo_num) - deallocate(a1) + !$OMP PARALLEL & + !$OMP DEFAULT(NONE) & + !$OMP PRIVATE(p, q, r) & + !$OMP SHARED(s, l, mo_num, mo_bi_ortho_tc_two_e_chemist, mo_r_coef, a_pqr) + !$OMP DO COLLAPSE(2) + do p = 1, mo_num + do q = 1, mo_num + do r = 1, mo_num + mo_bi_ortho_tc_two_e_chemist(p,q,r,s) = mo_bi_ortho_tc_two_e_chemist(p,q,r,s) + mo_r_coef(l,s) * a_pqr(p,q,r) + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL - call dgemm( 'T', 'N', mo_num*mo_num*mo_num, mo_num, ao_num, 1.d0 & - , a2(1,1,1,1), ao_num, mo_r_coef(1,1), ao_num & - , 0.d0, mo_bi_ortho_tc_two_e_chemist(1,1,1,1), mo_num*mo_num*mo_num) + enddo ! l + enddo ! s - deallocate(a2) + deallocate(a_jkp, a_kpq, a_pqr) + else + + print*, ' memory scale of TC ao -> mo: O(N4) ' + + allocate(a2(ao_num,ao_num,ao_num,mo_num)) + + call dgemm( 'T', 'N', ao_num*ao_num*ao_num, mo_num, ao_num, 1.d0 & + , ao_two_e_tc_tot(1,1,1,1), ao_num, mo_l_coef(1,1), ao_num & + , 0.d0, a2(1,1,1,1), ao_num*ao_num*ao_num) + + allocate(a1(ao_num,ao_num,mo_num,mo_num)) + + call dgemm( 'T', 'N', ao_num*ao_num*mo_num, mo_num, ao_num, 1.d0 & + , a2(1,1,1,1), ao_num, mo_r_coef(1,1), ao_num & + , 0.d0, a1(1,1,1,1), ao_num*ao_num*mo_num) + + deallocate(a2) + allocate(a2(ao_num,mo_num,mo_num,mo_num)) + + call dgemm( 'T', 'N', ao_num*mo_num*mo_num, mo_num, ao_num, 1.d0 & + , a1(1,1,1,1), ao_num, mo_l_coef(1,1), ao_num & + , 0.d0, a2(1,1,1,1), ao_num*mo_num*mo_num) + + deallocate(a1) + + call dgemm( 'T', 'N', mo_num*mo_num*mo_num, mo_num, ao_num, 1.d0 & + , a2(1,1,1,1), ao_num, mo_r_coef(1,1), ao_num & + , 0.d0, mo_bi_ortho_tc_two_e_chemist(1,1,1,1), mo_num*mo_num*mo_num) + + deallocate(a2) + + endif !allocate(a1(mo_num,ao_num,ao_num,ao_num)) !a1 = 0.d0 @@ -135,6 +192,10 @@ BEGIN_PROVIDER [double precision, mo_bi_ortho_tc_two_e_chemist, (mo_num, mo_num, !enddo !deallocate(a1) + call wall_time(t2) + print *, ' WALL TIME for PROVIDING mo_bi_ortho_tc_two_e_chemist (min)', (t2-t1)/60.d0 + call print_memory_usage() + END_PROVIDER ! --- diff --git a/plugins/local/non_h_ints_mu/total_tc_int.irp.f b/plugins/local/non_h_ints_mu/total_tc_int.irp.f index 9d3cf565..ba078d9b 100644 --- a/plugins/local/non_h_ints_mu/total_tc_int.irp.f +++ b/plugins/local/non_h_ints_mu/total_tc_int.irp.f @@ -201,6 +201,8 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n !$OMP END DO !$OMP END PARALLEL + call clear_ao_map() + if(tc_integ_type .eq. "numeric") then FREE int2_grad1_u12_ao_num int2_grad1_u12_square_ao_num endif diff --git a/src/tools/print_detweights.irp.f b/src/tools/print_detweights.irp.f new file mode 100644 index 00000000..d5b0f2c9 --- /dev/null +++ b/src/tools/print_detweights.irp.f @@ -0,0 +1,66 @@ +program print_detweights + + implicit none + + read_wf = .True. + touch read_wf + + call main() + +end + +! --- + +subroutine main() + + implicit none + integer :: i + integer :: degree + integer :: ios + integer, allocatable :: deg(:), ii(:), deg_sorted(:) + double precision, allocatable :: c(:) + + PROVIDE N_int + PROVIDE N_det + PROVIDE psi_det + PROVIDe psi_coef + + allocate(deg(N_det), ii(N_det), deg_sorted(N_det), c(N_det)) + + do i = 1, N_det + + call debug_det(psi_det(1,1,i), N_int) + call get_excitation_degree(psi_det(1,1,i), psi_det(1,1,1), degree, N_int) + + ii (i) = i + deg(i) = degree + c (i) = dabs(psi_coef(i,1)) + enddo + + call dsort(c, ii, N_det) + + do i = 1, N_det + deg_sorted(i) = deg(ii(i)) + enddo + + print *, ' saving psi' + + ! Writing output in binary format + open(unit=10, file="coef.bin", status="replace", action="write", iostat=ios, form="unformatted") + + if(ios /= 0) then + print *, ' Error opening file!' + stop + endif + + write(10) N_det + write(10) deg_sorted + write(10) c + + close(10) + + deallocate(deg, ii, deg_sorted, c) + +end + +