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