diff --git a/plugins/local/bi_ort_ints/semi_num_ints_mo.irp.f b/plugins/local/bi_ort_ints/semi_num_ints_mo.irp.f index 51f0cba4..77e4cb9b 100644 --- a/plugins/local/bi_ort_ints/semi_num_ints_mo.irp.f +++ b/plugins/local/bi_ort_ints/semi_num_ints_mo.irp.f @@ -107,8 +107,8 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_transp, (ao_num, ao_num, 3, integer :: i, j, ipoint double precision :: wall0, wall1 - print *, ' providing int2_grad1_u12_ao_transp ...' - call wall_time(wall0) + !print *, ' providing int2_grad1_u12_ao_transp ...' + !call wall_time(wall0) if(test_cycle_tc) then @@ -142,15 +142,15 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_transp, (ao_num, ao_num, 3, endif - call wall_time(wall1) - print *, ' wall time for int2_grad1_u12_ao_transp ', wall1 - wall0 - call print_memory_usage() + !call wall_time(wall1) + !print *, ' wall time for int2_grad1_u12_ao_transp (min) = ', (wall1 - wall0) / 60.d0 + !call print_memory_usage() END_PROVIDER ! --- -BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_transp, (mo_num, mo_num, 3, n_points_final_grid)] +BEGIN_PROVIDER [double precision, int2_grad1_u12_bimo_transp, (mo_num, mo_num, 3, n_points_final_grid)] implicit none integer :: ipoint @@ -159,7 +159,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_transp, (mo_num, mo_num, PROVIDE mo_l_coef mo_r_coef PROVIDE int2_grad1_u12_ao_transp - !print *, ' providing int2_grad1_u12_bimo_transp' + !print *, ' providing int2_grad1_u12_bimo_transp ...' !call wall_time(wall0) !$OMP PARALLEL & @@ -167,33 +167,35 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_transp, (mo_num, mo_num, !$OMP PRIVATE (ipoint) & !$OMP SHARED (n_points_final_grid,int2_grad1_u12_ao_transp,int2_grad1_u12_bimo_transp) !$OMP DO SCHEDULE (dynamic) - do ipoint = 1, n_points_final_grid - call ao_to_mo_bi_ortho( int2_grad1_u12_ao_transp (1,1,1,ipoint), size(int2_grad1_u12_ao_transp , 1) & - , int2_grad1_u12_bimo_transp(1,1,1,ipoint), size(int2_grad1_u12_bimo_transp, 1) ) - call ao_to_mo_bi_ortho( int2_grad1_u12_ao_transp (1,1,2,ipoint), size(int2_grad1_u12_ao_transp , 1) & - , int2_grad1_u12_bimo_transp(1,1,2,ipoint), size(int2_grad1_u12_bimo_transp, 1) ) - call ao_to_mo_bi_ortho( int2_grad1_u12_ao_transp (1,1,3,ipoint), size(int2_grad1_u12_ao_transp , 1) & - , int2_grad1_u12_bimo_transp(1,1,3,ipoint), size(int2_grad1_u12_bimo_transp, 1) ) - enddo + do ipoint = 1, n_points_final_grid + call ao_to_mo_bi_ortho( int2_grad1_u12_ao_transp (1,1,1,ipoint), size(int2_grad1_u12_ao_transp , 1) & + , int2_grad1_u12_bimo_transp(1,1,1,ipoint), size(int2_grad1_u12_bimo_transp, 1) ) + call ao_to_mo_bi_ortho( int2_grad1_u12_ao_transp (1,1,2,ipoint), size(int2_grad1_u12_ao_transp , 1) & + , int2_grad1_u12_bimo_transp(1,1,2,ipoint), size(int2_grad1_u12_bimo_transp, 1) ) + call ao_to_mo_bi_ortho( int2_grad1_u12_ao_transp (1,1,3,ipoint), size(int2_grad1_u12_ao_transp , 1) & + , int2_grad1_u12_bimo_transp(1,1,3,ipoint), size(int2_grad1_u12_bimo_transp, 1) ) + enddo !$OMP END DO !$OMP END PARALLEL + !FREE int2_grad1_u12_ao_transp + !call wall_time(wall1) - !print *, ' Wall time for providing int2_grad1_u12_bimo_transp',wall1 - wall0 + !print *, ' wall time for int2_grad1_u12_bimo_transp (min) =', (wall1 - wall0) / 60.d0 !call print_memory_usage() END_PROVIDER ! --- -BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_t, (n_points_final_grid, 3, mo_num, mo_num)] +BEGIN_PROVIDER [double precision, int2_grad1_u12_bimo_t, (n_points_final_grid, 3, mo_num, mo_num)] implicit none integer :: i, j, ipoint double precision :: wall0, wall1 !call wall_time(wall0) - !print *, ' Providing int2_grad1_u12_bimo_t ...' + !print *, ' providing int2_grad1_u12_bimo_t ...' PROVIDE mo_l_coef mo_r_coef PROVIDE int2_grad1_u12_bimo_transp @@ -211,17 +213,21 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_t, (n_points_final_grid, FREE int2_grad1_u12_bimo_transp !call wall_time(wall1) - !print *, ' wall time for int2_grad1_u12_bimo_t,', wall1 - wall0 + !print *, ' wall time for int2_grad1_u12_bimo_t (min) =', (wall1 - wall0) / 60.d0 !call print_memory_usage() END_PROVIDER ! --- -BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_t, (n_points_final_grid, 3, ao_num, ao_num)] +BEGIN_PROVIDER [double precision, int2_grad1_u12_ao_t, (n_points_final_grid, 3, ao_num, ao_num)] implicit none - integer :: i, j, ipoint + integer :: i, j, ipoint + double precision :: wall0, wall1 + + !call wall_time(wall0) + !print *, ' providing int2_grad1_u12_ao_t ...' PROVIDE int2_grad1_u12_ao @@ -235,6 +241,10 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_t, (n_points_final_grid, 3, enddo enddo + !call wall_time(wall1) + !print *, ' wall time for int2_grad1_u12_ao_t (min) =', (wall1 - wall0) / 60.d0 + !call print_memory_usage() + END_PROVIDER ! --- @@ -275,8 +285,8 @@ BEGIN_PROVIDER [ double precision, x_W_ki_bi_ortho_erf_rk, (n_points_final_grid, double precision :: xyz double precision :: wall0, wall1 - print*, ' providing x_W_ki_bi_ortho_erf_rk ...' - call wall_time(wall0) + !print*, ' providing x_W_ki_bi_ortho_erf_rk ...' + !call wall_time(wall0) !$OMP PARALLEL & !$OMP DEFAULT (NONE) & @@ -300,8 +310,8 @@ BEGIN_PROVIDER [ double precision, x_W_ki_bi_ortho_erf_rk, (n_points_final_grid, ! FREE mo_v_ki_bi_ortho_erf_rk_cst_mu_transp ! FREE mo_x_v_ki_bi_ortho_erf_rk_cst_mu_transp - call wall_time(wall1) - print *, ' time to provide x_W_ki_bi_ortho_erf_rk = ', wall1 - wall0 + !call wall_time(wall1) + !print *, ' time to provide x_W_ki_bi_ortho_erf_rk = ', wall1 - wall0 END_PROVIDER @@ -323,8 +333,8 @@ BEGIN_PROVIDER [ double precision, x_W_ki_bi_ortho_erf_rk_diag, (n_points_final_ double precision :: xyz double precision :: wall0, wall1 - print*,'providing x_W_ki_bi_ortho_erf_rk_diag ...' - call wall_time(wall0) + !print*,'providing x_W_ki_bi_ortho_erf_rk_diag ...' + !call wall_time(wall0) !$OMP PARALLEL & !$OMP DEFAULT (NONE) & @@ -343,8 +353,8 @@ BEGIN_PROVIDER [ double precision, x_W_ki_bi_ortho_erf_rk_diag, (n_points_final_ !$OMP END DO !$OMP END PARALLEL - call wall_time(wall1) - print*,'time to provide x_W_ki_bi_ortho_erf_rk_diag = ',wall1 - wall0 + !call wall_time(wall1) + !print*,'time to provide x_W_ki_bi_ortho_erf_rk_diag = ',wall1 - wall0 END_PROVIDER diff --git a/plugins/local/bi_ort_ints/three_body_ints_bi_ort.irp.f b/plugins/local/bi_ort_ints/three_body_ints_bi_ort.irp.f index 726e48ba..fd4a162f 100644 --- a/plugins/local/bi_ort_ints/three_body_ints_bi_ort.irp.f +++ b/plugins/local/bi_ort_ints/three_body_ints_bi_ort.irp.f @@ -168,7 +168,7 @@ subroutine give_integrals_3_body_bi_ort(n, l, k, m, j, i, integral) integral = integral + tmp * final_weight_at_r_vector(ipoint) enddo -end subroutine give_integrals_3_body_bi_ort +end ! --- 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 1e558038..e34a7b7b 100644 --- a/plugins/local/bi_ort_ints/total_twoe_pot.irp.f +++ b/plugins/local/bi_ort_ints/total_twoe_pot.irp.f @@ -16,10 +16,10 @@ double precision function bi_ortho_mo_ints(l, k, j, i) integer :: m, n, p, q bi_ortho_mo_ints = 0.d0 - do m = 1, ao_num - do p = 1, ao_num - do n = 1, ao_num - do q = 1, ao_num + do p = 1, ao_num + do m = 1, ao_num + do q = 1, ao_num + do n = 1, ao_num ! p1h1p2h2 l1 l2 r1 r2 bi_ortho_mo_ints += ao_two_e_tc_tot(n,q,m,p) * mo_l_coef(m,l) * mo_l_coef(n,k) * mo_r_coef(p,j) * mo_r_coef(q,i) enddo @@ -27,7 +27,7 @@ double precision function bi_ortho_mo_ints(l, k, j, i) enddo enddo -end function bi_ortho_mo_ints +end ! --- @@ -43,93 +43,68 @@ BEGIN_PROVIDER [double precision, mo_bi_ortho_tc_two_e_chemist, (mo_num, mo_num, integer :: i, j, k, l, m, n, p, q, s, r double precision :: t1, t2, tt1, tt2 double precision, allocatable :: a1(:,:,:,:), a2(:,:,:,:) - double precision, allocatable :: a_jkp(:,:,:), a_kpq(:,:,:), a_pqr(:,:,:) + double precision, allocatable :: a_jkp(:,:,:), a_kpq(:,:,:), ao_two_e_tc_tot_tmp(:,:,:) 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 if(ao_to_mo_tc_n3) then print*, ' memory scale of TC ao -> mo: O(N3) ' + if(.not.read_tc_integ) then + stop 'read_tc_integ needs to be set to true' + endif + 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)) + allocate(ao_two_e_tc_tot_tmp(ao_num,ao_num,ao_num)) + + open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/ao_two_e_tc_tot', action="read") call wall_time(tt1) - do s = 1, mo_num + mo_bi_ortho_tc_two_e_chemist(:,:,:,:) = 0.d0 + do l = 1, ao_num + read(11) ao_two_e_tc_tot_tmp(:,:,:) - mo_bi_ortho_tc_two_e_chemist(:,:,:,s) = 0.d0 - do l = 1, ao_num + do s = 1, 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 & + call dgemm( 'T', 'N', ao_num*ao_num, mo_num, ao_num, 1.d0 & + , ao_two_e_tc_tot_tmp(1,1,1), 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) - !$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, ao_num, mo_r_coef(l,s) & + , a_kpq(1,1,1), ao_num, mo_l_coef(1,1), ao_num & + , 1.d0, mo_bi_ortho_tc_two_e_chemist(1,1,1,s), mo_num*mo_num) - enddo ! l + enddo ! s - if(s == 2) then + if(l == 2) then call wall_time(tt2) print*, ' 1 / mo_num done in (min)', (tt2-tt1)/60.d0 print*, ' estimated time required (min)', dble(mo_num-1)*(tt2-tt1)/60.d0 - elseif(s == 11) then + elseif(l == 11) then call wall_time(tt2) print*, ' 10 / mo_num done in (min)', (tt2-tt1)/60.d0 print*, ' estimated time required (min)', dble(mo_num-10)*(tt2-tt1)/(60.d0*10.d0) - elseif(s == 26) then - call wall_time(tt2) - print*, ' 25 / mo_num done in (min)', (tt2-tt1)/60.d0 - print*, ' estimated time required (min)', dble(mo_num-25)*(tt2-tt1)/(60.d0*25.d0) - elseif(s == 51) then - call wall_time(tt2) - print*, ' 50 / mo_num done in (min)', (tt2-tt1)/60.d0 - print*, ' estimated time required (min)', dble(mo_num-50)*(tt2-tt1)/(60.d0*50.d0) - elseif(s == 101) then + elseif(l == 101) then call wall_time(tt2) print*, ' 100 / mo_num done in (min)', (tt2-tt1)/60.d0 print*, ' estimated time required (min)', dble(mo_num-100)*(tt2-tt1)/(60.d0*100.d0) - elseif(s == 201) then - call wall_time(tt2) - print*, ' 200 / mo_num done in (min)', (tt2-tt1)/60.d0 - print*, ' estimated time required (min)', dble(mo_num-200)*(tt2-tt1)/(60.d0*200.d0) - elseif(s == 501) then - call wall_time(tt2) - print*, ' 500 / mo_num done in (min)', (tt2-tt1)/60.d0 - print*, ' estimated time required (min)', dble(mo_num-500)*(tt2-tt1)/(60.d0*500.d0) endif + enddo ! l + close(11) - enddo ! s - - deallocate(a_jkp, a_kpq, a_pqr) + deallocate(a_jkp, a_kpq, ao_two_e_tc_tot_tmp) else @@ -141,6 +116,8 @@ BEGIN_PROVIDER [double precision, mo_bi_ortho_tc_two_e_chemist, (mo_num, mo_num, , 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) + FREE ao_two_e_tc_tot + 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 & diff --git a/plugins/local/non_h_ints_mu/tc_integ.irp.f b/plugins/local/non_h_ints_mu/tc_integ.irp.f index 775a9e4c..58e3db48 100644 --- a/plugins/local/non_h_ints_mu/tc_integ.irp.f +++ b/plugins/local/non_h_ints_mu/tc_integ.irp.f @@ -44,14 +44,92 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f elseif(tc_integ_type .eq. "numeric") then print *, ' Numerical integration over r1 and r2 will be performed' - - ! TODO combine 1shot & int2_grad1_u12_ao_num - PROVIDE int2_grad1_u12_ao_num - int2_grad1_u12_ao = int2_grad1_u12_ao_num + if(tc_save_mem) then - !PROVIDE int2_grad1_u12_ao_num_1shot - !int2_grad1_u12_ao = int2_grad1_u12_ao_num_1shot + integer :: n_blocks, n_rest, n_pass + integer :: i_blocks, i_rest, i_pass, ii + double precision :: mem, n_double + double precision, allocatable :: tmp(:,:,:), xx(:) + double precision, allocatable :: tmp_grad1_u12(:,:,:) + + PROVIDE final_weight_at_r_vector_extra aos_in_r_array_extra + + allocate(tmp(n_points_extra_final_grid,ao_num,ao_num), xx(n_points_extra_final_grid)) + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (j, i, jpoint) & + !$OMP SHARED (tmp, ao_num, n_points_extra_final_grid, final_weight_at_r_vector_extra, aos_in_r_array_extra_transp) + !$OMP DO COLLAPSE(2) + do j = 1, ao_num + do i = 1, ao_num + do jpoint = 1, n_points_extra_final_grid + tmp(jpoint,i,j) = final_weight_at_r_vector_extra(jpoint) * aos_in_r_array_extra_transp(jpoint,i) * aos_in_r_array_extra_transp(jpoint,j) + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + call total_memory(mem) + mem = max(1.d0, qp_max_mem - mem) + n_double = mem * 1.d8 + n_blocks = int(min(n_double / (n_points_extra_final_grid * 4.d0), 1.d0*n_points_final_grid)) + n_rest = int(mod(n_points_final_grid, n_blocks)) + n_pass = int((n_points_final_grid - n_rest) / n_blocks) + call write_int(6, n_pass, 'Number of passes') + call write_int(6, n_blocks, 'Size of the blocks') + call write_int(6, n_rest, 'Size of the last block') + allocate(tmp_grad1_u12(n_points_extra_final_grid,n_blocks,3)) + do i_pass = 1, n_pass + ii = (i_pass-1)*n_blocks + 1 + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i_blocks, ipoint) & + !$OMP SHARED (n_blocks, n_points_extra_final_grid, ii, final_grid_points, xx, tmp_grad1_u12) + !$OMP DO + do i_blocks = 1, n_blocks + ipoint = ii - 1 + i_blocks ! r1 + call get_grad1_u12_withsq_r1_seq(ipoint, n_points_extra_final_grid, tmp_grad1_u12(1,i_blocks,1), tmp_grad1_u12(1,i_blocks,2), tmp_grad1_u12(1,i_blocks,3), xx(1)) + enddo + !$OMP END DO + !$OMP END PARALLEL + do m = 1, 3 + call dgemm( "T", "N", ao_num*ao_num, n_blocks, n_points_extra_final_grid, 1.d0 & + , tmp(1,1,1), n_points_extra_final_grid, tmp_grad1_u12(1,1,m), n_points_extra_final_grid & + , 0.d0, int2_grad1_u12_ao(1,1,ii,m), ao_num*ao_num) + enddo + enddo + deallocate(tmp_grad1_u12) + if(n_rest .gt. 0) then + allocate(tmp_grad1_u12(n_points_extra_final_grid,n_rest,3)) + ii = n_pass*n_blocks + 1 + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i_rest, ipoint) & + !$OMP SHARED (n_rest, n_points_extra_final_grid, ii, final_grid_points, xx, tmp_grad1_u12) + !$OMP DO + do i_rest = 1, n_rest + ipoint = ii - 1 + i_rest ! r1 + call get_grad1_u12_withsq_r1_seq(ipoint, n_points_extra_final_grid, tmp_grad1_u12(1,i_rest,1), tmp_grad1_u12(1,i_rest,2), tmp_grad1_u12(1,i_rest,3), xx(1)) + enddo + !$OMP END DO + !$OMP END PARALLEL + do m = 1, 3 + call dgemm( "T", "N", ao_num*ao_num, n_rest, n_points_extra_final_grid, 1.d0 & + , tmp(1,1,1), n_points_extra_final_grid, tmp_grad1_u12(1,1,m), n_points_extra_final_grid & + , 0.d0, int2_grad1_u12_ao(1,1,ii,m), ao_num*ao_num) + enddo + deallocate(tmp_grad1_u12) + endif + deallocate(tmp,xx) + + else + ! TODO combine 1shot & int2_grad1_u12_ao_num + PROVIDE int2_grad1_u12_ao_num + int2_grad1_u12_ao = int2_grad1_u12_ao_num + !PROVIDE int2_grad1_u12_ao_num_1shot + !int2_grad1_u12_ao = int2_grad1_u12_ao_num_1shot + endif elseif(tc_integ_type .eq. "semi-analytic") then @@ -177,13 +255,88 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_p print *, ' Numerical integration over r1 and r2 will be performed' - ! TODO combine 1shot & int2_grad1_u12_square_ao_num + if(tc_save_mem) then - PROVIDE int2_grad1_u12_square_ao_num - int2_grad1_u12_square_ao = int2_grad1_u12_square_ao_num + integer :: n_blocks, n_rest, n_pass + integer :: i_blocks, i_rest, i_pass, ii + double precision :: mem, n_double + double precision, allocatable :: tmp(:,:,:), xx(:,:,:) + double precision, allocatable :: tmp_grad1_u12_squared(:,:) - !PROVIDE int2_grad1_u12_square_ao_num_1shot - !int2_grad1_u12_square_ao = int2_grad1_u12_square_ao_num_1shot + PROVIDE final_weight_at_r_vector_extra aos_in_r_array_extra + + allocate(tmp(n_points_extra_final_grid,ao_num,ao_num)) + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (j, i, jpoint) & + !$OMP SHARED (tmp, ao_num, n_points_extra_final_grid, final_weight_at_r_vector_extra, aos_in_r_array_extra_transp) + !$OMP DO COLLAPSE(2) + do j = 1, ao_num + do i = 1, ao_num + do jpoint = 1, n_points_extra_final_grid + tmp(jpoint,i,j) = final_weight_at_r_vector_extra(jpoint) * aos_in_r_array_extra_transp(jpoint,i) * aos_in_r_array_extra_transp(jpoint,j) + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + call total_memory(mem) + mem = max(1.d0, qp_max_mem - mem) + n_double = mem * 1.d8 + n_blocks = int(min(n_double / (n_points_extra_final_grid * 4.d0), 1.d0*n_points_final_grid)) + n_rest = int(mod(n_points_final_grid, n_blocks)) + n_pass = int((n_points_final_grid - n_rest) / n_blocks) + call write_int(6, n_pass, 'Number of passes') + call write_int(6, n_blocks, 'Size of the blocks') + call write_int(6, n_rest, 'Size of the last block') + allocate(tmp_grad1_u12_squared(n_points_extra_final_grid,n_blocks), xx(n_points_extra_final_grid,n_blocks,3)) + do i_pass = 1, n_pass + ii = (i_pass-1)*n_blocks + 1 + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i_blocks, ipoint) & + !$OMP SHARED (n_blocks, n_points_extra_final_grid, ii, xx, final_grid_points, tmp_grad1_u12_squared) + !$OMP DO + do i_blocks = 1, n_blocks + ipoint = ii - 1 + i_blocks ! r1 + call get_grad1_u12_withsq_r1_seq(ipoint, n_points_extra_final_grid, xx(1,i_blocks,1), xx(1,i_blocks,2), xx(1,i_blocks,3), tmp_grad1_u12_squared(1,i_blocks)) + enddo + !$OMP END DO + !$OMP END PARALLEL + call dgemm( "T", "N", ao_num*ao_num, n_blocks, n_points_extra_final_grid, -0.5d0 & + , tmp(1,1,1), n_points_extra_final_grid, tmp_grad1_u12_squared(1,1), n_points_extra_final_grid & + , 0.d0, int2_grad1_u12_square_ao(1,1,ii), ao_num*ao_num) + enddo + deallocate(tmp_grad1_u12_squared, xx) + if(n_rest .gt. 0) then + ii = n_pass*n_blocks + 1 + allocate(tmp_grad1_u12_squared(n_points_extra_final_grid,n_rest), xx(n_points_extra_final_grid,n_rest,3)) + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i_rest, ipoint) & + !$OMP SHARED (n_rest, n_points_extra_final_grid, ii, xx, final_grid_points, tmp_grad1_u12_squared) + !$OMP DO + do i_rest = 1, n_rest + ipoint = ii - 1 + i_rest ! r1 + call get_grad1_u12_withsq_r1_seq(ipoint, n_points_extra_final_grid, xx(1,i_rest,1), xx(1,i_rest,2), xx(1,i_rest,3), tmp_grad1_u12_squared(1,i_rest)) + enddo + !$OMP END DO + !$OMP END PARALLEL + call dgemm( "T", "N", ao_num*ao_num, n_rest, n_points_extra_final_grid, -0.5d0 & + , tmp(1,1,1), n_points_extra_final_grid, tmp_grad1_u12_squared(1,1), n_points_extra_final_grid & + , 0.d0, int2_grad1_u12_square_ao(1,1,ii), ao_num*ao_num) + deallocate(tmp_grad1_u12_squared, xx) + endif + deallocate(tmp) + + else + + ! TODO combine 1shot & int2_grad1_u12_square_ao_num + PROVIDE int2_grad1_u12_square_ao_num + int2_grad1_u12_square_ao = int2_grad1_u12_square_ao_num + !PROVIDE int2_grad1_u12_square_ao_num_1shot + !int2_grad1_u12_square_ao = int2_grad1_u12_square_ao_num_1shot + endif elseif(tc_integ_type .eq. "semi-analytic") then 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 b8379006..a1bbd6e0 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 @@ -55,7 +55,9 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n print*, ' Reading ao_two_e_tc_tot from ', trim(ezfio_filename) // '/work/ao_two_e_tc_tot' open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/ao_two_e_tc_tot', action="read") - read(11) ao_two_e_tc_tot + do i = 1, ao_num + read(11) ao_two_e_tc_tot(:,:,:,i) + enddo close(11) else @@ -67,7 +69,7 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n PROVIDE int2_grad1_u12_square_ao - if(tc_save_mem) then + if(tc_save_mem_loops) then print*, ' LOOPS are used to evaluate Hermitian part of ao_two_e_tc_tot ...' @@ -176,7 +178,7 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n PROVIDE int2_grad1_u12_ao - if(tc_save_mem) then + if(tc_save_mem_loops) then print*, ' LOOPS are used to evaluate non-Hermitian part of ao_two_e_tc_tot ...' @@ -241,7 +243,6 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n deallocate(c_mat) end if - !FREE int2_grad1_u12_ao if(tc_integ_type .eq. "semi-analytic") then FREE int2_grad1_u2e_ao @@ -264,48 +265,52 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n print*, ' adding ERI to ao_two_e_tc_tot ...' - !$OMP PARALLEL DEFAULT(NONE) & - !$OMP PRIVATE(i, j, k, l, integ_zero, integ_val) & - !$OMP SHARED(ao_num, ao_two_e_tc_tot) - !$OMP DO COLLAPSE(4) - do j = 1, ao_num - do l = 1, ao_num - do i = 1, ao_num - do k = 1, ao_num - integ_zero = ao_two_e_integral_zero(i,j,k,l) - if(.not. integ_zero) then - ! i,k : r1 j,l : r2 - integ_val = ao_two_e_integral(i,k,j,l) - ao_two_e_tc_tot(k,i,l,j) = ao_two_e_tc_tot(k,i,l,j) + integ_val - endif + if(tc_save_mem) then + print*, ' ao_integrals_map will not be used' + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(i, j, k, l, integ_zero, integ_val) & + !$OMP SHARED(ao_num, ao_two_e_tc_tot) + !$OMP DO COLLAPSE(4) + do j = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do k = 1, ao_num + integ_zero = ao_two_e_integral_zero(i,j,k,l) + if(.not. integ_zero) then + ! i,k : r1 j,l : r2 + integ_val = ao_two_e_integral(i,k,j,l) + ao_two_e_tc_tot(k,i,l,j) = ao_two_e_tc_tot(k,i,l,j) + integ_val + endif + enddo enddo enddo enddo - enddo - !$OMP END DO - !$OMP END PARALLEL + !$OMP END DO + !$OMP END PARALLEL + else + print*, ' ao_integrals_map will be used' + PROVIDE ao_integrals_map + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP SHARED(ao_num, ao_two_e_tc_tot, ao_integrals_map) & + !$OMP PRIVATE(i, j, k, l) + !$OMP DO COLLAPSE(4) + do j = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do k = 1, ao_num + ! < 1:i, 2:j | 1:k, 2:l > + ao_two_e_tc_tot(k,i,l,j) = ao_two_e_tc_tot(k,i,l,j) + get_ao_two_e_integral(i, j, k, l, ao_integrals_map) + enddo + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + !call clear_ao_map() + FREE ao_integrals_map + endif - !PROVIDE ao_integrals_map - !!$OMP PARALLEL DEFAULT(NONE) & - !!$OMP SHARED(ao_num, ao_two_e_tc_tot, ao_integrals_map) & - !!$OMP PRIVATE(i, j, k, l) - !!$OMP DO COLLAPSE(4) - !do j = 1, ao_num - ! do l = 1, ao_num - ! do i = 1, ao_num - ! do k = 1, ao_num - ! ! < 1:i, 2:j | 1:k, 2:l > - ! ao_two_e_tc_tot(k,i,l,j) = ao_two_e_tc_tot(k,i,l,j) + get_ao_two_e_integral(i, j, k, l, ao_integrals_map) - ! enddo - ! enddo - ! enddo - !enddo - !!$OMP END DO - !!$OMP END PARALLEL - !!call clear_ao_map() - !FREE ao_integrals_map - - if(tc_integ_type .eq. "numeric") then + if((tc_integ_type .eq. "numeric") .and. (.not. tc_save_mem)) then FREE int2_grad1_u12_ao_num int2_grad1_u12_square_ao_num endif @@ -315,7 +320,9 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n print*, ' Saving ao_two_e_tc_tot in ', trim(ezfio_filename) // '/work/ao_two_e_tc_tot' open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/ao_two_e_tc_tot', action="write") call ezfio_set_work_empty(.False.) - write(11) ao_two_e_tc_tot + do i = 1, ao_num + write(11) ao_two_e_tc_tot(:,:,:,i) + enddo close(11) call ezfio_set_tc_keywords_io_tc_integ('Read') endif diff --git a/plugins/local/tc_keywords/EZFIO.cfg b/plugins/local/tc_keywords/EZFIO.cfg index 24362796..70169189 100644 --- a/plugins/local/tc_keywords/EZFIO.cfg +++ b/plugins/local/tc_keywords/EZFIO.cfg @@ -286,10 +286,16 @@ doc: If |true|, memory scale of TC ao -> mo: O(N3) interface: ezfio,provider,ocaml default: False -[tc_save_mem] +[tc_save_mem_loops] type: logical doc: If |true|, use loops to save memory TC interface: ezfio,provider,ocaml default: False +[tc_save_mem] +type: logical +doc: If |true|, more calc but less mem +interface: ezfio,provider,ocaml +default: False + diff --git a/plugins/local/tc_scf/fock_3e_bi_ortho_cs.irp.f b/plugins/local/tc_scf/fock_3e_bi_ortho_cs.irp.f index 0b883865..8fd5e5b6 100644 --- a/plugins/local/tc_scf/fock_3e_bi_ortho_cs.irp.f +++ b/plugins/local/tc_scf/fock_3e_bi_ortho_cs.irp.f @@ -9,7 +9,7 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_cs, (mo_num, mo_num)] double precision :: loc_1, loc_2, loc_3 double precision, allocatable :: Okappa(:), Jkappa(:,:) double precision, allocatable :: tmp_omp_d1(:), tmp_omp_d2(:,:) - double precision, allocatable :: tmp_1(:,:), tmp_2(:,:,:,:) + double precision, allocatable :: tmp_1(:,:), tmp_2(:,:,:,:), tmp_22(:,:,:) double precision, allocatable :: tmp_3(:,:,:), tmp_4(:,:,:) PROVIDE mo_l_coef mo_r_coef @@ -63,17 +63,13 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_cs, (mo_num, mo_num)] allocate(tmp_1(n_points_final_grid,4)) do ipoint = 1, n_points_final_grid - loc_1 = 2.d0 * Okappa(ipoint) - tmp_1(ipoint,1) = loc_1 * Jkappa(ipoint,1) tmp_1(ipoint,2) = loc_1 * Jkappa(ipoint,2) tmp_1(ipoint,3) = loc_1 * Jkappa(ipoint,3) - tmp_1(ipoint,4) = Okappa(ipoint) enddo - !$OMP PARALLEL & !$OMP DEFAULT (NONE) & !$OMP PRIVATE (ipoint, i, j, loc_1, tmp_omp_d2) & @@ -112,58 +108,81 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_cs, (mo_num, mo_num)] ! --- - allocate(tmp_2(n_points_final_grid,4,mo_num,mo_num)) + if(tc_save_mem) then - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (ipoint, a, b) & - !$OMP SHARED (n_points_final_grid, mo_num, & - !$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, & - !$OMP int2_grad1_u12_bimo_t, final_weight_at_r_vector, & - !$OMP tmp_2) - !$OMP DO COLLAPSE(2) - do a = 1, mo_num - do b = 1, mo_num - do ipoint = 1, n_points_final_grid - tmp_2(ipoint,1,b,a) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,b,a) - tmp_2(ipoint,2,b,a) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,b,a) - tmp_2(ipoint,3,b,a) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,b,a) - enddo - enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (ipoint, a, b, i) & - !$OMP SHARED (n_points_final_grid, mo_num, elec_beta_num, & - !$OMP final_weight_at_r_vector, int2_grad1_u12_bimo_t, & - !$OMP tmp_2) - !$OMP DO COLLAPSE(2) - do a = 1, mo_num - do b = 1, mo_num - tmp_2(:,4,b,a) = 0.d0 - do i = 1, elec_beta_num + allocate(tmp_22(n_points_final_grid,4,mo_num)) + do a = 1, mo_num + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, b, i) & + !$OMP SHARED (n_points_final_grid, mo_num, elec_beta_num, a, & + !$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, & + !$OMP int2_grad1_u12_bimo_t, final_weight_at_r_vector, & + !$OMP tmp_22) + !$OMP DO + do b = 1, mo_num do ipoint = 1, n_points_final_grid - tmp_2(ipoint,4,b,a) -= final_weight_at_r_vector(ipoint) * ( int2_grad1_u12_bimo_t(ipoint,1,b,i) * int2_grad1_u12_bimo_t(ipoint,1,i,a) & - + int2_grad1_u12_bimo_t(ipoint,2,b,i) * int2_grad1_u12_bimo_t(ipoint,2,i,a) & - + int2_grad1_u12_bimo_t(ipoint,3,b,i) * int2_grad1_u12_bimo_t(ipoint,3,i,a) ) + tmp_22(ipoint,1,b) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,b,a) + tmp_22(ipoint,2,b) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,b,a) + tmp_22(ipoint,3,b) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,b,a) + enddo + tmp_22(:,4,b) = 0.d0 + do i = 1, elec_beta_num + do ipoint = 1, n_points_final_grid + tmp_22(ipoint,4,b) -= final_weight_at_r_vector(ipoint) * ( int2_grad1_u12_bimo_t(ipoint,1,b,i) * int2_grad1_u12_bimo_t(ipoint,1,i,a) & + + int2_grad1_u12_bimo_t(ipoint,2,b,i) * int2_grad1_u12_bimo_t(ipoint,2,i,a) & + + int2_grad1_u12_bimo_t(ipoint,3,b,i) * int2_grad1_u12_bimo_t(ipoint,3,i,a) ) + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + call dgemv( 'T', 4*n_points_final_grid, mo_num, -2.d0 & + , tmp_22(1,1,1), size(tmp_22, 1) * size(tmp_22, 2) & + , tmp_1(1,1), 1 & + , 0.d0, fock_3e_uhf_mo_cs(1,a), 1) + enddo + deallocate(tmp_22) + + else + + allocate(tmp_2(n_points_final_grid,4,mo_num,mo_num)) + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, a, b, i) & + !$OMP SHARED (n_points_final_grid, mo_num, elec_beta_num, & + !$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, & + !$OMP int2_grad1_u12_bimo_t, final_weight_at_r_vector, & + !$OMP tmp_2) + !$OMP DO COLLAPSE(2) + do a = 1, mo_num + do b = 1, mo_num + do ipoint = 1, n_points_final_grid + tmp_2(ipoint,1,b,a) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,b,a) + tmp_2(ipoint,2,b,a) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,b,a) + tmp_2(ipoint,3,b,a) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,b,a) + enddo + tmp_2(:,4,b,a) = 0.d0 + do i = 1, elec_beta_num + do ipoint = 1, n_points_final_grid + tmp_2(ipoint,4,b,a) -= final_weight_at_r_vector(ipoint) * ( int2_grad1_u12_bimo_t(ipoint,1,b,i) * int2_grad1_u12_bimo_t(ipoint,1,i,a) & + + int2_grad1_u12_bimo_t(ipoint,2,b,i) * int2_grad1_u12_bimo_t(ipoint,2,i,a) & + + int2_grad1_u12_bimo_t(ipoint,3,b,i) * int2_grad1_u12_bimo_t(ipoint,3,i,a) ) + enddo enddo enddo enddo - enddo - !$OMP END DO - !$OMP END PARALLEL + !$OMP END DO + !$OMP END PARALLEL + call dgemv( 'T', 4*n_points_final_grid, mo_num*mo_num, -2.d0 & + , tmp_2(1,1,1,1), size(tmp_2, 1) * size(tmp_2, 2) & + , tmp_1(1,1), 1 & + , 0.d0, fock_3e_uhf_mo_cs(1,1), 1) + deallocate(tmp_2) - ! --- + endif - call dgemv( 'T', 4*n_points_final_grid, mo_num*mo_num, -2.d0 & - , tmp_2(1,1,1,1), size(tmp_2, 1) * size(tmp_2, 2) & - , tmp_1(1,1), 1 & - , 0.d0, fock_3e_uhf_mo_cs(1,1), 1) - - deallocate(tmp_1, tmp_2) + deallocate(tmp_1) ! --- @@ -272,7 +291,7 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_cs, (mo_num, mo_num)] ! --- !call wall_time(tf) - !print *, ' total Wall time for fock_3e_uhf_mo_cs =', tf - ti + !print *, ' total Wall time for fock_3e_uhf_mo_cs =', (tf - ti) / 60.d0 END_PROVIDER diff --git a/plugins/local/tc_scf/fock_3e_bi_ortho_uhf.irp.f b/plugins/local/tc_scf/fock_3e_bi_ortho_uhf.irp.f index 63a1e162..47ee5b48 100644 --- a/plugins/local/tc_scf/fock_3e_bi_ortho_uhf.irp.f +++ b/plugins/local/tc_scf/fock_3e_bi_ortho_uhf.irp.f @@ -32,7 +32,7 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_a, (mo_num, mo_num)] endif !call wall_time(tf) - !print *, ' Wall time for fock_3e_uhf_mo_a =', tf - ti + !print *, ' Wall time for fock_3e_uhf_mo_a (min) =', (tf - ti)/60.d0 END_PROVIDER diff --git a/plugins/local/tc_scf/fock_tc.irp.f b/plugins/local/tc_scf/fock_tc.irp.f index 282f9873..d3ddb8ad 100644 --- a/plugins/local/tc_scf/fock_tc.irp.f +++ b/plugins/local/tc_scf/fock_tc.irp.f @@ -175,7 +175,7 @@ END_PROVIDER ! --- -BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_alpha, (mo_num, mo_num) ] +BEGIN_PROVIDER [double precision, Fock_matrix_tc_mo_alpha, (mo_num, mo_num)] BEGIN_DOC ! Total alpha TC Fock matrix : h_c + Two-e^TC terms on the MO basis diff --git a/plugins/local/tc_scf/write_ao_2e_tc_integ.irp.f b/plugins/local/tc_scf/write_ao_2e_tc_integ.irp.f new file mode 100644 index 00000000..7ce57578 --- /dev/null +++ b/plugins/local/tc_scf/write_ao_2e_tc_integ.irp.f @@ -0,0 +1,58 @@ +! --- + +program write_ao_2e_tc_integ + + implicit none + + PROVIDE j1e_type + PROVIDE j2e_type + + print *, ' j1e_type = ', j1e_type + print *, ' j2e_type = ', j2e_type + + my_grid_becke = .True. + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + + call write_int(6, my_n_pt_r_grid, 'radial external grid over') + call write_int(6, my_n_pt_a_grid, 'angular external grid over') + + if(tc_integ_type .eq. "numeric") then + my_extra_grid_becke = .True. + PROVIDE tc_grid2_a tc_grid2_r + my_n_pt_r_extra_grid = tc_grid2_r + my_n_pt_a_extra_grid = tc_grid2_a + touch my_extra_grid_becke my_n_pt_r_extra_grid my_n_pt_a_extra_grid + + call write_int(6, my_n_pt_r_extra_grid, 'radial internal grid over') + call write_int(6, my_n_pt_a_extra_grid, 'angular internal grid over') + endif + + call main() + +end + +! --- + +subroutine main() + + implicit none + + PROVIDE io_tc_integ + + print*, 'io_tc_integ = ', io_tc_integ + + if(io_tc_integ .ne. "Write") then + print*, 'io_tc_integ != Write' + print*, io_tc_integ + stop + endif + + PROVIDE ao_two_e_tc_tot + +end + +! --- + diff --git a/src/mol_properties/multi_s_dipole_moment.irp.f b/src/mol_properties/multi_s_dipole_moment.irp.f index f21e08cd..c7216a61 100644 --- a/src/mol_properties/multi_s_dipole_moment.irp.f +++ b/src/mol_properties/multi_s_dipole_moment.irp.f @@ -102,12 +102,28 @@ END_PROVIDER &BEGIN_PROVIDER [double precision, multi_s_z_dipole_moment_eigenval, (N_states)] implicit none + double precision, allocatable :: eigval(:), eigvec(:,:), A(:,:) PROVIDE multi_s_x_dipole_moment multi_s_y_dipole_moment multi_s_z_dipole_moment - call lapack_diag(multi_s_x_dipole_moment_eigenval(1), multi_s_x_dipole_moment_eigenvec(1,1), multi_s_x_dipole_moment(1,1), N_states, N_states) - call lapack_diag(multi_s_y_dipole_moment_eigenval(1), multi_s_y_dipole_moment_eigenvec(1,1), multi_s_y_dipole_moment(1,1), N_states, N_states) - call lapack_diag(multi_s_z_dipole_moment_eigenval(1), multi_s_z_dipole_moment_eigenvec(1,1), multi_s_z_dipole_moment(1,1), N_states, N_states) + allocate(A(N_states,N_states), eigvec(N_states,N_states), eigval(N_states)) + + A = multi_s_x_dipole_moment + call lapack_diag(eigval(1), eigvec(1,1), A(1,1), N_states, N_states) + multi_s_x_dipole_moment_eigenval = eigval + multi_s_x_dipole_moment_eigenvec = eigvec + + A = multi_s_y_dipole_moment + call lapack_diag(eigval(1), eigvec(1,1), A(1,1), N_states, N_states) + multi_s_y_dipole_moment_eigenval = eigval + multi_s_y_dipole_moment_eigenvec = eigvec + + A = multi_s_z_dipole_moment + call lapack_diag(eigval(1), eigvec(1,1), A(1,1), N_states, N_states) + multi_s_z_dipole_moment_eigenval = eigval + multi_s_z_dipole_moment_eigenvec = eigvec + + deallocate(A, eigvec, eigval) END_PROVIDER