From 7039831efcc521119e0c4fb221507eb270bedf1c Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 22 Feb 2021 16:25:05 +0100 Subject: [PATCH 1/5] Added Abdallah's file --- devel/svdwf/Evar_TruncSVD.irp.f | 167 ++++++++++++++++++++++++++++++++ devel/svdwf/NEED | 1 + devel/svdwf/random_svd.irp.f | 14 ++- 3 files changed, 177 insertions(+), 5 deletions(-) create mode 100644 devel/svdwf/Evar_TruncSVD.irp.f diff --git a/devel/svdwf/Evar_TruncSVD.irp.f b/devel/svdwf/Evar_TruncSVD.irp.f new file mode 100644 index 0000000..4ae793d --- /dev/null +++ b/devel/svdwf/Evar_TruncSVD.irp.f @@ -0,0 +1,167 @@ +program Evar_TruncSVD + implicit none + BEGIN_DOC +! study energy variation with truncated SVD + END_DOC + read_wf = .True. + TOUCH read_wf + ! !!! + call run() + ! !!! +end + + + + +subroutine run + implicit none + include 'constants.include.F' + double precision, allocatable :: A(:,:), U(:,:), V(:,:), D(:) + integer :: r, i, j, k, l, m, n, iter, iter_max + double precision, allocatable :: Z(:,:), P(:,:), Yt(:,:), UYt(:,:), r1(:,:) + ! !!! + m = n_det_alpha_unique + n = n_det_beta_unique + r = n + print *, 'matrix:', m,'x',n + print *, 'N det:', N_det + print *, 'rank = ', r + iter_max = 20 + ! !!! + allocate( Z(m,r) , P(n,r) ) ! Z(m,r) = A(m,n) @ P(n,r) + Z(:,:) = 0.d0 + ! !!! + ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! + ! first we apply a RSVD for a pre-fixed rank (r) + ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! + !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l,r1) + allocate(r1(N_det,2)) + !$OMP DO + do l=1,r + call random_number(r1) + r1(:,1) = dsqrt(-2.d0*dlog(r1(:,1))) + r1(:,1) = r1(:,1) * dcos(dtwo_pi*r1(:,2)) + 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) + enddo + enddo + !$OMP END DO + deallocate(r1) + !$OMP END PARALLEL + ! !!! + ! Power iterations + ! !!! + do iter=1,iter_max + ! !!! + print *, 'Power iteration ', iter, '/', 20 + ! !!! + ! P(n,r) = At(n,m) @ Z(m,r) + ! !!! + !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l) + !$OMP DO + do l=1,r + P(:,l) = 0.d0 + do k=1,N_det + i = psi_bilinear_matrix_rows(k) + j = psi_bilinear_matrix_columns(k) + P(j,l) = P(j,l) + psi_bilinear_matrix_values(k,1) * Z(i,l) + enddo + enddo + !$OMP END DO + ! !!! + ! Z(m,r) = A(m,n) @ P(n,r) + ! !!! + !$OMP BARRIER + !$OMP DO + do l=1,r + Z(:,l) = 0.d0 + 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) * P(j,l) + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + ! !!! + ! Compute QR: at return: Q is Z(m,r) + ! !!! + call ortho_qr(Z,size(Z,1),m,r) + ! !!! + enddo + ! !!! + ! Y(r,n) = Zt(r,m) @ A(m,n) or Yt(n,r) = At(n,m) @ Z(m,r) + ! !!! + allocate(Yt(n,r)) + !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l) + !$OMP DO + do l=1,r + do k=1,n + Yt(k,l) = 0.d0 + enddo + do k=1,N_det + i = psi_bilinear_matrix_rows(k) + j = psi_bilinear_matrix_columns(k) + Yt(j,l) = Yt(j,l) + Z(i,l) * psi_bilinear_matrix_values(k,1) + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + ! !!! + ! Y = UY @ D @ Vt or Yt = V @ Dt @ UYt + ! !!! + allocate(D(r),V(n,r), UYt(r,r)) + ! !!! + 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) or U = Z @ (UYt).T + ! !!! + 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) + ! !!! + !do i=1,r + ! print *, i, real(D(i)), real(D(i)**2), real(sum(D(1:i)**2)) + ! if (D(i) < 1.d-15) then + ! k = i + ! exit + ! endif + !enddo + !print *, 'threshold: ', 2.858 * D(k/2) + ! !!! + ! Build the new determinant: U @ D @ Vt + ! !!! + !!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l) + !!$OMP DO + !! + !print *, 'ok 1' + !N_det = m * n + !print *, 'ok 11' + !TOUCH N_det + !psi_bilinear_matrix_values(:,1) = 0.d0 + !TOUCH psi_bilinear_matrix_values + ! print *, size(psi_bilinear_matrix_values,1), size(D), size(U,1), size(U,2), size(V,1), size(V,2) + print*, PSI_energy(1) + nuclear_repulsion + psi_bilinear_matrix(:,:,:) = 0.d0 + do r = 1, n + call generate_all_alpha_beta_det_products + do i = 1, N_det_beta_unique + do j = 1, N_det_alpha_unique + psi_bilinear_matrix(j,i,1) = 0.d0 + do l = 1, r + psi_bilinear_matrix(j,i,1) = psi_bilinear_matrix(j,i,1) + D(l) * U(j,l) * V(i,l) + enddo + enddo + enddo + TOUCH psi_bilinear_matrix + call update_wf_of_psi_bilinear_matrix(.False.) + print*, r, PSI_energy(1) + nuclear_repulsion !CI_energy(1) + enddo + !!$OMP END DO + !!$OMP END PARALLEL + deallocate(U,D,V) + ! !!! +end diff --git a/devel/svdwf/NEED b/devel/svdwf/NEED index d3d4d2c..8d89a45 100644 --- a/devel/svdwf/NEED +++ b/devel/svdwf/NEED @@ -1 +1,2 @@ determinants +davidson_undressed diff --git a/devel/svdwf/random_svd.irp.f b/devel/svdwf/random_svd.irp.f index 9ac9119..2ae76dc 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 @@ -98,7 +102,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) From 0f890c872d639bf1448bc3a64335f8fd2618c66f Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 8 Mar 2021 12:54:34 +0100 Subject: [PATCH 2/5] add hmo.qp to import_integrals --- .../import_integrals_mo.irp.f | 149 +++++++++++------- devel/svdwf/svdwf.irp.f | 4 + stable/utilities/truncate_wf.irp.f | 8 +- 3 files changed, 97 insertions(+), 64 deletions(-) diff --git a/devel/import_integrals/import_integrals_mo.irp.f b/devel/import_integrals/import_integrals_mo.irp.f index 966fe87..7784c87 100644 --- a/devel/import_integrals/import_integrals_mo.irp.f +++ b/devel/import_integrals/import_integrals_mo.irp.f @@ -25,6 +25,8 @@ subroutine run ! ! Vmo.qp : electron-nucleus potential ! +! hmo.qp : core hamiltonian. If hmo exists, the other 1-e files are not read +! ! Wmo.qp : electron repulsion integrals ! END_DOC @@ -53,68 +55,95 @@ subroutine run call ezfio_set_nuclei_io_nuclear_repulsion('Read') endif - A = 0.d0 - i = 1 - j = 1 - iunit = getunitandopen('Tmo.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 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') + logical :: exists + inquire(file='hmo.qp',exist=exists) + if (exists) then + A = 0.d0 + i = 1 + j = 1 + iunit = getunitandopen('hmo.qp','r') + do while(.true.) + read (iunit,*,end=23) i,j, integral + if (i<0 .or. i>mo_num) then + print *, i + stop 'i out of bounds in hmo.qp' + endif + if (j<0 .or. j>mo_num) then + print *, j + stop 'j out of bounds in hmo.qp' + endif + A(i,j) = integral + enddo + 23 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') - 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') + else + A = 0.d0 + i = 1 + j = 1 + iunit = getunitandopen('Tmo.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 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') - 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('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') + + 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') + + endif iunit = getunitandopen('Wmo.qp','r') n_integrals=0 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) From 2399a3ecaa30ff836ac8f9504e349856b7882eed Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 12 Mar 2021 13:13:49 +0100 Subject: [PATCH 3/5] Added Hmo.qp --- .../export_integrals_mo.irp.f | 1 + .../import_integrals_mo.irp.f | 176 +++++++++++------- devel/svdwf/Evar_TruncSVD.irp.f | 5 +- 3 files changed, 110 insertions(+), 72 deletions(-) diff --git a/devel/import_integrals/export_integrals_mo.irp.f b/devel/import_integrals/export_integrals_mo.irp.f index 292a2bc..20f2ac9 100644 --- a/devel/import_integrals/export_integrals_mo.irp.f +++ b/devel/import_integrals/export_integrals_mo.irp.f @@ -61,6 +61,7 @@ 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 diff --git a/devel/import_integrals/import_integrals_mo.irp.f b/devel/import_integrals/import_integrals_mo.irp.f index 966fe87..0d94bec 100644 --- a/devel/import_integrals/import_integrals_mo.irp.f +++ b/devel/import_integrals/import_integrals_mo.irp.f @@ -19,29 +19,33 @@ subroutine run ! ! Tmo.qp : kinetic energy integrals ! -! Smo.qp : overlap matrix -! ! Pmo.qp : pseudopotential integrals ! ! Vmo.qp : electron-nucleus potential + +! Hmo.qp : one-electron integrals. If Hmo is present, don't try to read V,P,T ! ! Wmo.qp : electron repulsion integrals ! END_DOC - integer :: iunit - integer :: getunitandopen + integer :: iunit + integer :: getunitandopen - integer ::i,j,k,l - double precision :: integral - double precision, allocatable :: A(:,:) + integer :: i,j,k,l + double precision :: integral + double precision, allocatable :: A(:,:) - integer :: n_integrals - integer(key_kind), allocatable :: buffer_i(:) + integer :: n_integrals + logical :: exists + integer(key_kind), allocatable :: buffer_i(:) real(integral_kind), allocatable :: buffer_values(:) allocate(buffer_i(mo_num**3), buffer_values(mo_num**3)) allocate (A(mo_num,mo_num)) + PROVIDE mo_two_e_integrals_in_map + + ! Nuclear repulsion A(1,1) = huge(1.d0) iunit = getunitandopen('E.qp','r') @@ -53,68 +57,102 @@ subroutine run call ezfio_set_nuclei_io_nuclear_repulsion('Read') endif - A = 0.d0 - i = 1 - j = 1 - iunit = getunitandopen('Tmo.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 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') - 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') + ! One-electron integrals - 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') + exists = .False. + inquire(file='Hmo.qp',exist=exists) + if (exists) then + + A = 0.d0 + i = 1 + j = 1 + iunit = getunitandopen('Hmo.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 Hmo.qp' + endif + if (j<0 .or. j>mo_num) then + print *, j + stop 'j out of bounds in Hmo.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') + + else + + A = 0.d0 + i = 1 + j = 1 + iunit = getunitandopen('Tmo.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 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') + + 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') + + 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') + + end if + + ! Two-electron integrals iunit = getunitandopen('Wmo.qp','r') n_integrals=0 diff --git a/devel/svdwf/Evar_TruncSVD.irp.f b/devel/svdwf/Evar_TruncSVD.irp.f index 4ae793d..3476d0d 100644 --- a/devel/svdwf/Evar_TruncSVD.irp.f +++ b/devel/svdwf/Evar_TruncSVD.irp.f @@ -158,10 +158,9 @@ subroutine run enddo TOUCH psi_bilinear_matrix call update_wf_of_psi_bilinear_matrix(.False.) - print*, r, PSI_energy(1) + nuclear_repulsion !CI_energy(1) + print*, r, PSI_energy(1) + nuclear_repulsion, s2_values(1) !CI_energy(1) + call save_wavefunction() enddo - !!$OMP END DO - !!$OMP END PARALLEL deallocate(U,D,V) ! !!! end From 04443346847b75751c329a48a58d8601e993d364 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 17 Mar 2021 14:00:56 +0100 Subject: [PATCH 4/5] Bug fix --- devel/import_integrals/export_integrals_mo.irp.f | 2 ++ 1 file changed, 2 insertions(+) diff --git a/devel/import_integrals/export_integrals_mo.irp.f b/devel/import_integrals/export_integrals_mo.irp.f index 292a2bc..5e47145 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 @@ -55,6 +56,7 @@ subroutine run double precision, external :: get_two_e_integral + allocate (A(mo_num,mo_num)) call ezfio_set_nuclei_nuclear_repulsion(nuclear_repulsion) From ec6c19596cddb6cb075e53fc495e9a213089e1c2 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 17 Mar 2021 14:08:24 +0100 Subject: [PATCH 5/5] Bug fix in import_integrals_mo --- .../export_integrals_mo.irp.f | 9 +-- .../import_integrals_mo.irp.f | 63 ++++++++++--------- 2 files changed, 40 insertions(+), 32 deletions(-) diff --git a/devel/import_integrals/export_integrals_mo.irp.f b/devel/import_integrals/export_integrals_mo.irp.f index ba8a387..be0d558 100644 --- a/devel/import_integrals/export_integrals_mo.irp.f +++ b/devel/import_integrals/export_integrals_mo.irp.f @@ -54,7 +54,7 @@ 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)) @@ -68,9 +68,10 @@ subroutine run 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 Hmo.qp' + 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 Hmo.qp' + stop 'j out of bounds in H_mo.qp' endif A(i,j) = integral enddo @@ -83,22 +84,27 @@ subroutine run 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') else + call ezfio_set_mo_one_e_ints_io_mo_one_e_integrals('None') + A = 0.d0 i = 1 j = 1 - iunit = getunitandopen('Tmo.qp','r') + 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 Tmo.qp' + 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 Tmo.qp' + stop 'j out of bounds in T_mo.qp' endif A(i,j) = integral enddo @@ -110,16 +116,16 @@ subroutine run A = 0.d0 i = 1 j = 1 - iunit = getunitandopen('Pmo.qp','r') + 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 Pmo.qp' + 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 Pmo.qp' + stop 'j out of bounds in P_mo.qp' endif A(i,j) = integral enddo @@ -131,16 +137,16 @@ subroutine run A = 0.d0 i = 1 j = 1 - iunit = getunitandopen('Vmo.qp','r') + 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 Vmo.qp' + 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 Vmo.qp' + stop 'j out of bounds in V_mo.qp' endif A(i,j) = integral enddo @@ -154,7 +160,7 @@ subroutine run ! Two-electron integrals - iunit = getunitandopen('Wmo.qp','r') + iunit = getunitandopen('W_mo.qp','r') n_integrals=0 i = 1 j = 1 @@ -165,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 @@ -192,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 + +