1
0
mirror of https://gitlab.com/scemama/qp_plugins_scemama.git synced 2025-04-19 06:50:11 +02:00

add hmo.qp to import_integrals

This commit is contained in:
Anthony Scemama 2021-03-08 12:54:34 +01:00
parent 8371b5318f
commit 0f890c872d
3 changed files with 97 additions and 64 deletions

View File

@ -25,6 +25,8 @@ subroutine run
! !
! Vmo.qp : electron-nucleus potential ! 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 ! Wmo.qp : electron repulsion integrals
! !
END_DOC END_DOC
@ -53,68 +55,95 @@ subroutine run
call ezfio_set_nuclei_io_nuclear_repulsion('Read') call ezfio_set_nuclei_io_nuclear_repulsion('Read')
endif endif
A = 0.d0 logical :: exists
i = 1 inquire(file='hmo.qp',exist=exists)
j = 1 if (exists) then
iunit = getunitandopen('Tmo.qp','r') A = 0.d0
do i = 1
read (iunit,*,end=10) i,j, integral j = 1
if (i<0 .or. i>mo_num) then iunit = getunitandopen('hmo.qp','r')
print *, i do while(.true.)
stop 'i out of bounds in Tmo.qp' read (iunit,*,end=23) i,j, integral
endif if (i<0 .or. i>mo_num) then
if (j<0 .or. j>mo_num) then print *, i
print *, j stop 'i out of bounds in hmo.qp'
stop 'j out of bounds in Tmo.qp' endif
endif if (j<0 .or. j>mo_num) then
A(i,j) = integral print *, j
enddo stop 'j out of bounds in hmo.qp'
10 continue endif
close(iunit) A(i,j) = integral
call ezfio_set_mo_one_e_ints_mo_integrals_kinetic(A) enddo
call ezfio_set_mo_one_e_ints_io_mo_integrals_kinetic('Read') 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 else
i = 1 A = 0.d0
j = 1 i = 1
iunit = getunitandopen('Pmo.qp','r') j = 1
do iunit = getunitandopen('Tmo.qp','r')
read (iunit,*,end=14) i,j, integral do
if (i<0 .or. i>mo_num) then read (iunit,*,end=10) i,j, integral
print *, i if (i<0 .or. i>mo_num) then
stop 'i out of bounds in Pmo.qp' print *, i
endif stop 'i out of bounds in Tmo.qp'
if (j<0 .or. j>mo_num) then endif
print *, j if (j<0 .or. j>mo_num) then
stop 'j out of bounds in Pmo.qp' print *, j
endif stop 'j out of bounds in Tmo.qp'
A(i,j) = integral endif
enddo A(i,j) = integral
14 continue enddo
close(iunit) 10 continue
call ezfio_set_mo_one_e_ints_mo_integrals_pseudo(A) close(iunit)
call ezfio_set_mo_one_e_ints_io_mo_integrals_pseudo('Read') 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 A = 0.d0
i = 1 i = 1
j = 1 j = 1
iunit = getunitandopen('Vmo.qp','r') iunit = getunitandopen('Pmo.qp','r')
do do
read (iunit,*,end=12) i,j, integral read (iunit,*,end=14) i,j, integral
if (i<0 .or. i>mo_num) then if (i<0 .or. i>mo_num) then
print *, i print *, i
stop 'i out of bounds in Vmo.qp' stop 'i out of bounds in Pmo.qp'
endif endif
if (j<0 .or. j>mo_num) then if (j<0 .or. j>mo_num) then
print *, j print *, j
stop 'j out of bounds in Vmo.qp' stop 'j out of bounds in Pmo.qp'
endif endif
A(i,j) = integral A(i,j) = integral
enddo enddo
12 continue 14 continue
close(iunit) close(iunit)
call ezfio_set_mo_one_e_ints_mo_integrals_n_e(A) call ezfio_set_mo_one_e_ints_mo_integrals_pseudo(A)
call ezfio_set_mo_one_e_ints_io_mo_integrals_n_e('Read') 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') iunit = getunitandopen('Wmo.qp','r')
n_integrals=0 n_integrals=0

View File

@ -53,6 +53,10 @@ subroutine run
print *, 'threshold: ', 2.858 * D(k/2) print *, 'threshold: ', 2.858 * D(k/2)
print *, 'Entropy : ', entropy 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 ! do i=1,n_det_alpha_unique
! print '(I6,4(X,F12.8))', i, U(i,1:4) ! print '(I6,4(X,F12.8))', i, U(i,1:4)
! enddo ! enddo

View File

@ -55,16 +55,16 @@ subroutine routine_s2
integer :: i,j,k integer :: i,j,k
double precision :: accu(N_states) double precision :: accu(N_states)
print *, 'Weights of the SOP' print *, 'Weights of the CFG'
do i=1,N_det 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 enddo
print*, 'Min weight of the occupation pattern ?' print*, 'Min weight of the occupation pattern ?'
read(5,*) wmin read(5,*) wmin
ndet_max = 0 ndet_max = 0
do i=1,N_det 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 ndet_max = ndet_max+1
enddo enddo
@ -73,7 +73,7 @@ subroutine routine_s2
accu = 0.d0 accu = 0.d0
k=0 k=0
do i = 1, N_det 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 k = k+1
do j = 1, N_int do j = 1, N_int
psi_det_tmp(j,1,k) = psi_det(j,1,i) psi_det_tmp(j,1,k) = psi_det(j,1,i)