diff --git a/devel/import_integrals/export_integrals_mo.irp.f b/devel/import_integrals/export_integrals_mo.irp.f index 292a2bc..be0d558 100644 --- a/devel/import_integrals/export_integrals_mo.irp.f +++ b/devel/import_integrals/export_integrals_mo.irp.f @@ -1,4 +1,5 @@ program export_integrals_mo + PROVIDE mo_two_e_integrals_in_map call run end @@ -53,7 +54,8 @@ subroutine run integer :: n_integrals integer(key_kind) :: idx - double precision, external :: get_two_e_integral + double precision, external :: mo_two_e_integral + allocate (A(mo_num,mo_num)) call ezfio_set_nuclei_nuclear_repulsion(nuclear_repulsion) @@ -61,13 +63,15 @@ subroutine run call write_2d('T_mo.qp',mo_kinetic_integrals) call write_2d('P_mo.qp',mo_pseudo_integrals) call write_2d('V_mo.qp',mo_integrals_n_e) + call write_2d('H_mo.qp',mo_one_e_integrals) iunit = getunitandopen('W_mo.qp','w') do l=1,mo_num do k=1,mo_num - do j=1,mo_num - do i=1,mo_num - integral = get_two_e_integral(i,j,k,l,mo_integrals_map) + do j=l,mo_num + do i=k,mo_num + if (i==j .and. kmo_num) then - print *, i - stop 'i out of bounds in Tmo.qp' - endif - if (j<0 .or. j>mo_num) then - print *, j - stop 'j out of bounds in Tmo.qp' - endif - A(i,j) = integral - enddo - 10 continue - close(iunit) - call ezfio_set_mo_one_e_ints_mo_integrals_kinetic(A) - call ezfio_set_mo_one_e_ints_io_mo_integrals_kinetic('Read') + ! One-electron integrals - A = 0.d0 - i = 1 - j = 1 - iunit = getunitandopen('Pmo.qp','r') - do - read (iunit,*,end=14) i,j, integral - if (i<0 .or. i>mo_num) then - print *, i - stop 'i out of bounds in Pmo.qp' - endif - if (j<0 .or. j>mo_num) then - print *, j - stop 'j out of bounds in Pmo.qp' - endif - A(i,j) = integral - enddo - 14 continue - close(iunit) - call ezfio_set_mo_one_e_ints_mo_integrals_pseudo(A) - call ezfio_set_mo_one_e_ints_io_mo_integrals_pseudo('Read') + exists = .False. + inquire(file='H_mo.qp',exist=exists) + if (exists) then - A = 0.d0 - i = 1 - j = 1 - iunit = getunitandopen('Vmo.qp','r') - do - read (iunit,*,end=12) i,j, integral - if (i<0 .or. i>mo_num) then - print *, i - stop 'i out of bounds in Vmo.qp' - endif - if (j<0 .or. j>mo_num) then - print *, j - stop 'j out of bounds in Vmo.qp' - endif - A(i,j) = integral - enddo - 12 continue - close(iunit) - call ezfio_set_mo_one_e_ints_mo_integrals_n_e(A) - call ezfio_set_mo_one_e_ints_io_mo_integrals_n_e('Read') + A = 0.d0 + i = 1 + j = 1 + iunit = getunitandopen('H_mo.qp','r') + do + read (iunit,*,end=8) i,j, integral + if (i<0 .or. i>mo_num) then + print *, i + stop 'i out of bounds in H_mo.qp' + endif + if (j<0 .or. j>mo_num) then + print *, j + stop 'j out of bounds in H_mo.qp' + endif + A(i,j) = integral + enddo + 8 continue + close(iunit) + call ezfio_set_mo_one_e_ints_mo_one_e_integrals(A) + call ezfio_set_mo_one_e_ints_io_mo_one_e_integrals('Read') + call ezfio_set_mo_one_e_ints_io_mo_integrals_kinetic('None') + call ezfio_set_mo_one_e_ints_io_mo_integrals_pseudo('None') + call ezfio_set_mo_one_e_ints_io_mo_integrals_n_e('None') - iunit = getunitandopen('Wmo.qp','r') + else + + call ezfio_set_mo_one_e_ints_io_mo_one_e_integrals('None') + + A = 0.d0 + i = 1 + j = 1 + iunit = getunitandopen('T_mo.qp','r') + do + read (iunit,*,end=10) i,j, integral + if (i<0 .or. i>mo_num) then + print *, i + stop 'i out of bounds in T_mo.qp' + endif + if (j<0 .or. j>mo_num) then + print *, j + stop 'j out of bounds in T_mo.qp' + endif + A(i,j) = integral + enddo + 10 continue + close(iunit) + call ezfio_set_mo_one_e_ints_mo_integrals_kinetic(A) + call ezfio_set_mo_one_e_ints_io_mo_integrals_kinetic('Read') + + A = 0.d0 + i = 1 + j = 1 + iunit = getunitandopen('P_mo.qp','r') + do + read (iunit,*,end=14) i,j, integral + if (i<0 .or. i>mo_num) then + print *, i + stop 'i out of bounds in P_mo.qp' + endif + if (j<0 .or. j>mo_num) then + print *, j + stop 'j out of bounds in P_mo.qp' + endif + A(i,j) = integral + enddo + 14 continue + close(iunit) + call ezfio_set_mo_one_e_ints_mo_integrals_pseudo(A) + call ezfio_set_mo_one_e_ints_io_mo_integrals_pseudo('Read') + + A = 0.d0 + i = 1 + j = 1 + iunit = getunitandopen('V_mo.qp','r') + do + read (iunit,*,end=12) i,j, integral + if (i<0 .or. i>mo_num) then + print *, i + stop 'i out of bounds in V_mo.qp' + endif + if (j<0 .or. j>mo_num) then + print *, j + stop 'j out of bounds in V_mo.qp' + endif + A(i,j) = integral + enddo + 12 continue + close(iunit) + call ezfio_set_mo_one_e_ints_mo_integrals_n_e(A) + call ezfio_set_mo_one_e_ints_io_mo_integrals_n_e('Read') + + end if + + + ! Two-electron integrals + + iunit = getunitandopen('W_mo.qp','r') n_integrals=0 i = 1 j = 1 @@ -127,26 +171,25 @@ subroutine run read (iunit,*,end=13) i,j,k,l, integral if (i<0 .or. i>mo_num) then print *, i - stop 'i out of bounds in Wmo.qp' + stop 'i out of bounds in W_mo.qp' endif if (j<0 .or. j>mo_num) then print *, j - stop 'j out of bounds in Wmo.qp' + stop 'j out of bounds in W_mo.qp' endif if (k<0 .or. k>mo_num) then print *, k - stop 'k out of bounds in Wmo.qp' + stop 'k out of bounds in W_mo.qp' endif if (l<0 .or. l>mo_num) then print *, l - stop 'l out of bounds in Wmo.qp' + stop 'l out of bounds in W_mo.qp' endif n_integrals += 1 call mo_two_e_integrals_index(i, j, k, l, buffer_i(n_integrals) ) buffer_values(n_integrals) = integral if (n_integrals == size(buffer_i)) then - call insert_into_mo_integrals_map(n_integrals, buffer_i, buffer_values, & - real(mo_integrals_threshold,integral_kind)) + call map_append(mo_integrals_map, buffer_i, buffer_values, n_integrals) n_integrals = 0 endif enddo @@ -154,13 +197,15 @@ subroutine run close(iunit) if (n_integrals > 0) then - call insert_into_mo_integrals_map(n_integrals, buffer_i, buffer_values, & - real(mo_integrals_threshold,integral_kind)) + call map_append(mo_integrals_map, buffer_i, buffer_values, n_integrals) endif + call map_sort(mo_integrals_map) call map_unique(mo_integrals_map) call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints',mo_integrals_map) call ezfio_set_mo_two_e_ints_io_mo_two_e_integrals('Read') end + + diff --git a/devel/svdwf/Evar_TruncSVD.irp.f b/devel/svdwf/Evar_TruncSVD.irp.f index 1a72dc5..aa480e1 100644 --- a/devel/svdwf/Evar_TruncSVD.irp.f +++ b/devel/svdwf/Evar_TruncSVD.irp.f @@ -44,14 +44,14 @@ subroutine run print *, 'matrix:', m,'x',n print *, 'N det:', N_det ! !!! - r_init = 78 !!!42 - delta_r = 1 + r_init = 78 + delta_r = 5 ! !!! PowerIt_max = 10 nb_oversamp = 10 ! !!! It_max = 10 - error_thr = 1.d-3 !not stable if less thant 0.001 + error_thr = 1.d-3 !! don't touche it but rather increase r_init ! !!! ! !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! ~ !!! ! ! !!! ~ Rank Revealing Randomized SVD ~ !!! ! diff --git a/devel/svdwf/random_svd.irp.f b/devel/svdwf/random_svd.irp.f index 9bfb50b..5819519 100644 --- a/devel/svdwf/random_svd.irp.f +++ b/devel/svdwf/random_svd.irp.f @@ -22,13 +22,17 @@ subroutine run allocate(Z(m,r)) - ! Z(m,r) = A(m,n).P(n,r) - Z(:,:) = 0.d0 + ! Z(m,r) = A(m,n).P(n,r) + do j=1,r + do i=1,m + Z(i,j) = 0.d0 + enddo + enddo allocate(P(n,r)) !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l,r1) allocate(r1(N_det,2)) - !$OMP DO + !$OMP DO do l=1,r call random_number(r1) r1(:,1) = dsqrt(-2.d0*dlog(r1(:,1))) @@ -36,7 +40,7 @@ subroutine run do k=1,N_det i = psi_bilinear_matrix_rows(k) j = psi_bilinear_matrix_columns(k) - Z(i,l) = Z(i,l) + psi_bilinear_matrix_values(k,1) * r1(k,1) + Z(i,l) = Z(i,l) + psi_bilinear_matrix_values(k,1) * r1(k,1) enddo enddo !$OMP END DO @@ -100,7 +104,7 @@ subroutine run call svd(Yt,size(Yt,1),V,size(V,1),D,UYt,size(UYt,1),n,r) deallocate(Yt) - ! U(m,r) = Z(m,r).UY(r,r) + ! U(m,r) = Z(m,r).UY(r,r) allocate(U(m,r)) call dgemm('N','T',m,r,r,1.d0,Z,size(Z,1),UYt,size(UYt,1),0.d0,U,size(U,1)) deallocate(UYt,Z) diff --git a/devel/svdwf/svdwf.irp.f b/devel/svdwf/svdwf.irp.f index 1f67022..c694b3e 100644 --- a/devel/svdwf/svdwf.irp.f +++ b/devel/svdwf/svdwf.irp.f @@ -53,6 +53,10 @@ subroutine run print *, 'threshold: ', 2.858 * D(k/2) print *, 'Entropy : ', entropy + call ezfio_set_spindeterminants_psi_svd_alpha(U) + call ezfio_set_spindeterminants_psi_svd_beta (Vt) + call ezfio_set_spindeterminants_psi_svd_coefs(D) + ! do i=1,n_det_alpha_unique ! print '(I6,4(X,F12.8))', i, U(i,1:4) ! enddo diff --git a/stable/utilities/truncate_wf.irp.f b/stable/utilities/truncate_wf.irp.f index 6fe00ef..f1904d5 100644 --- a/stable/utilities/truncate_wf.irp.f +++ b/stable/utilities/truncate_wf.irp.f @@ -55,16 +55,16 @@ subroutine routine_s2 integer :: i,j,k double precision :: accu(N_states) - print *, 'Weights of the SOP' + print *, 'Weights of the CFG' do i=1,N_det - print *, i, real(weight_occ_pattern(det_to_occ_pattern(i),:)), real(sum(weight_occ_pattern(det_to_occ_pattern(i),:))) + print *, i, real(weight_configuration(det_to_configuration(i),:)), real(sum(weight_configuration(det_to_configuration(i),:))) enddo print*, 'Min weight of the occupation pattern ?' read(5,*) wmin ndet_max = 0 do i=1,N_det - if (maxval(weight_occ_pattern( det_to_occ_pattern(i),:)) < wmin) cycle + if (maxval(weight_configuration( det_to_configuration(i),:)) < wmin) cycle ndet_max = ndet_max+1 enddo @@ -73,7 +73,7 @@ subroutine routine_s2 accu = 0.d0 k=0 do i = 1, N_det - if (maxval(weight_occ_pattern( det_to_occ_pattern(i),:)) < wmin) cycle + if (maxval(weight_configuration( det_to_configuration(i),:)) < wmin) cycle k = k+1 do j = 1, N_int psi_det_tmp(j,1,k) = psi_det(j,1,i)