From 0f890c872d639bf1448bc3a64335f8fd2618c66f Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 8 Mar 2021 12:54:34 +0100 Subject: [PATCH] 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)