1
0
mirror of https://gitlab.com/scemama/qp_plugins_scemama.git synced 2024-11-09 07:33:41 +01:00

Compare commits

..

No commits in common. "ec6c19596cddb6cb075e53fc495e9a213089e1c2" and "c3890acae73a13b4460f1dc998c6c8db984778f7" have entirely different histories.

2 changed files with 32 additions and 42 deletions

View File

@ -1,5 +1,4 @@
program export_integrals_mo
PROVIDE mo_two_e_integrals_in_map
call run
end
@ -54,8 +53,7 @@ subroutine run
integer :: n_integrals
integer(key_kind) :: idx
double precision, external :: mo_two_e_integral
double precision, external :: get_two_e_integral
allocate (A(mo_num,mo_num))
call ezfio_set_nuclei_nuclear_repulsion(nuclear_repulsion)
@ -68,10 +66,9 @@ subroutine run
iunit = getunitandopen('W_mo.qp','w')
do l=1,mo_num
do k=1,mo_num
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)
do j=1,mo_num
do i=1,mo_num
integral = get_two_e_integral(i,j,k,l,mo_integrals_map)
if (integral /= 0.d0) then
write (iunit,'(4(I5,2X),E22.15)') i,j,k,l, integral
endif

View File

@ -17,15 +17,15 @@ subroutine run
!
! E.qp : Contains the nuclear repulsion energy
!
! T_mo.qp : kinetic energy integrals
! Tmo.qp : kinetic energy integrals
!
! P_mo.qp : pseudopotential integrals
! Pmo.qp : pseudopotential integrals
!
! V_mo.qp : electron-nucleus potential
! Vmo.qp : electron-nucleus potential
!
! H_mo.qp : core hamiltonian. If hmo exists, the other 1-e files are not read
! Hmo.qp : core hamiltonian. If hmo exists, the other 1-e files are not read
!
! W_mo.qp : electron repulsion integrals
! Wmo.qp : electron repulsion integrals
!
END_DOC
@ -43,8 +43,7 @@ subroutine run
allocate(buffer_i(mo_num**3), buffer_values(mo_num**3))
allocate (A(mo_num,mo_num))
PROVIDE mo_integrals_map
PROVIDE mo_integrals_threshold
PROVIDE mo_two_e_integrals_in_map
! Nuclear repulsion
@ -61,22 +60,22 @@ subroutine run
! One-electron integrals
exists = .False.
inquire(file='H_mo.qp',exist=exists)
inquire(file='Hmo.qp',exist=exists)
if (exists) then
A = 0.d0
i = 1
j = 1
iunit = getunitandopen('H_mo.qp','r')
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 H_mo.qp'
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 H_mo.qp'
stop 'j out of bounds in Hmo.qp'
endif
A(i,j) = integral
enddo
@ -84,27 +83,22 @@ 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('T_mo.qp','r')
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 T_mo.qp'
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 T_mo.qp'
stop 'j out of bounds in Tmo.qp'
endif
A(i,j) = integral
enddo
@ -116,16 +110,16 @@ subroutine run
A = 0.d0
i = 1
j = 1
iunit = getunitandopen('P_mo.qp','r')
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 P_mo.qp'
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 P_mo.qp'
stop 'j out of bounds in Pmo.qp'
endif
A(i,j) = integral
enddo
@ -137,16 +131,16 @@ subroutine run
A = 0.d0
i = 1
j = 1
iunit = getunitandopen('V_mo.qp','r')
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 V_mo.qp'
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 V_mo.qp'
stop 'j out of bounds in Vmo.qp'
endif
A(i,j) = integral
enddo
@ -160,7 +154,7 @@ subroutine run
! Two-electron integrals
iunit = getunitandopen('W_mo.qp','r')
iunit = getunitandopen('Wmo.qp','r')
n_integrals=0
i = 1
j = 1
@ -171,25 +165,26 @@ 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 W_mo.qp'
stop 'i out of bounds in Wmo.qp'
endif
if (j<0 .or. j>mo_num) then
print *, j
stop 'j out of bounds in W_mo.qp'
stop 'j out of bounds in Wmo.qp'
endif
if (k<0 .or. k>mo_num) then
print *, k
stop 'k out of bounds in W_mo.qp'
stop 'k out of bounds in Wmo.qp'
endif
if (l<0 .or. l>mo_num) then
print *, l
stop 'l out of bounds in W_mo.qp'
stop 'l out of bounds in Wmo.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 map_append(mo_integrals_map, buffer_i, buffer_values, n_integrals)
call insert_into_mo_integrals_map(n_integrals, buffer_i, buffer_values, &
real(mo_integrals_threshold,integral_kind))
n_integrals = 0
endif
enddo
@ -197,15 +192,13 @@ subroutine run
close(iunit)
if (n_integrals > 0) then
call map_append(mo_integrals_map, buffer_i, buffer_values, n_integrals)
call insert_into_mo_integrals_map(n_integrals, buffer_i, buffer_values, &
real(mo_integrals_threshold,integral_kind))
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