scemama
/
qp_plugins_scemama
mirror of https://gitlab.com/scemama/qp_plugins_scemama.git
1
0
Fork 0
Code Issues Releases Wiki Activity

Bug fix in import_integrals_mo

This commit is contained in:
Anthony Scemama 2021-03-17 14:08:24 +01:00
parent 564ac4b089
commit ec6c19596c
2 changed files with 40 additions and 32 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
devel/import_integrals/export_integrals_mo.irp.f
View File

@ -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. k<l) cycle
integral = mo_two_e_integral(i,j,k,l)
if (integral /= 0.d0) then
write (iunit,'(4(I5,2X),E22.15)') i,j,k,l, integral
endif

63
devel/import_integrals/import_integrals_mo.irp.f
View File

@ -17,15 +17,15 @@ subroutine run
!
! E.qp : Contains the nuclear repulsion energy
!
! Tmo.qp : kinetic energy integrals
! T_mo.qp : kinetic energy integrals
!
! Pmo.qp : pseudopotential integrals
! P_mo.qp : pseudopotential integrals
!
! Vmo.qp : electron-nucleus potential
! V_mo.qp : electron-nucleus potential
!
! Hmo.qp : core hamiltonian. If hmo exists, the other 1-e files are not read
! H_mo.qp : core hamiltonian. If hmo exists, the other 1-e files are not read
!
! Wmo.qp : electron repulsion integrals
! W_mo.qp : electron repulsion integrals
!
END_DOC
@ -43,7 +43,8 @@ subroutine run
allocate(buffer_i(mo_num**3), buffer_values(mo_num**3))
allocate (A(mo_num,mo_num))
PROVIDE mo_two_e_integrals_in_map
PROVIDE mo_integrals_map
PROVIDE mo_integrals_threshold
! Nuclear repulsion
@ -60,22 +61,22 @@ subroutine run
! One-electron integrals
exists = .False.
inquire(file='Hmo.qp',exist=exists)
inquire(file='H_mo.qp',exist=exists)
if (exists) then
A = 0.d0
i = 1
j = 1
iunit = getunitandopen('Hmo.qp','r')
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 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