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

Compare commits

base: scemama:ec6c19596cddb6cb075e53fc495e9a213089e1c2
Branches Tags
scemama:master
..
compare: scemama:c3890acae73a13b4460f1dc998c6c8db984778f7
Branches Tags
scemama:master

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

2 changed files with 32 additions and 42 deletions
Show Stats Expand all files Collapse all files

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

@ -1,5 +1,4 @@
program export_integrals_mo program export_integrals_mo
PROVIDE mo_two_e_integrals_in_map
call run call run
end end
@ -54,8 +53,7 @@ subroutine run
integer :: n_integrals integer :: n_integrals
integer(key_kind) :: idx integer(key_kind) :: idx
double precision, external :: mo_two_e_integral double precision, external :: get_two_e_integral
allocate (A(mo_num,mo_num)) allocate (A(mo_num,mo_num))
call ezfio_set_nuclei_nuclear_repulsion(nuclear_repulsion) call ezfio_set_nuclei_nuclear_repulsion(nuclear_repulsion)
@ -68,10 +66,9 @@ subroutine run
iunit = getunitandopen('W_mo.qp','w') iunit = getunitandopen('W_mo.qp','w')
do l=1,mo_num do l=1,mo_num
do k=1,mo_num do k=1,mo_num
do j=l,mo_num do j=1,mo_num
do i=k,mo_num do i=1,mo_num
if (i==j .and. k<l) cycle integral = get_two_e_integral(i,j,k,l,mo_integrals_map)
integral = mo_two_e_integral(i,j,k,l)
if (integral /= 0.d0) then if (integral /= 0.d0) then
write (iunit,'(4(I5,2X),E22.15)') i,j,k,l, integral write (iunit,'(4(I5,2X),E22.15)') i,j,k,l, integral
endif 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 ! 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 END_DOC
@ -43,8 +43,7 @@ subroutine run
allocate(buffer_i(mo_num**3), buffer_values(mo_num**3)) allocate(buffer_i(mo_num**3), buffer_values(mo_num**3))
allocate (A(mo_num,mo_num)) allocate (A(mo_num,mo_num))
PROVIDE mo_integrals_map PROVIDE mo_two_e_integrals_in_map
PROVIDE mo_integrals_threshold
! Nuclear repulsion ! Nuclear repulsion
@ -61,22 +60,22 @@ subroutine run
! One-electron integrals ! One-electron integrals
exists = .False. exists = .False.
inquire(file='H_mo.qp',exist=exists) inquire(file='Hmo.qp',exist=exists)
if (exists) then if (exists) then
A = 0.d0 A = 0.d0
i = 1 i = 1
j = 1 j = 1
iunit = getunitandopen('H_mo.qp','r') iunit = getunitandopen('Hmo.qp','r')
do do
read (iunit,*,end=8) i,j, integral read (iunit,*,end=8) 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 H_mo.qp' stop 'i out of bounds in Hmo.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 H_mo.qp' stop 'j out of bounds in Hmo.qp'
endif endif
A(i,j) = integral A(i,j) = integral
enddo enddo
@ -84,27 +83,22 @@ subroutine run
close(iunit) close(iunit)
call ezfio_set_mo_one_e_ints_mo_one_e_integrals(A) 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_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 else
call ezfio_set_mo_one_e_ints_io_mo_one_e_integrals('None')
A = 0.d0 A = 0.d0
i = 1 i = 1
j = 1 j = 1
iunit = getunitandopen('T_mo.qp','r') iunit = getunitandopen('Tmo.qp','r')
do do
read (iunit,*,end=10) i,j, integral read (iunit,*,end=10) 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 T_mo.qp' stop 'i out of bounds in Tmo.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 T_mo.qp' stop 'j out of bounds in Tmo.qp'
endif endif
A(i,j) = integral A(i,j) = integral
enddo enddo
@ -116,16 +110,16 @@ subroutine run
A = 0.d0 A = 0.d0
i = 1 i = 1
j = 1 j = 1
iunit = getunitandopen('P_mo.qp','r') iunit = getunitandopen('Pmo.qp','r')
do do
read (iunit,*,end=14) 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 P_mo.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 P_mo.qp' stop 'j out of bounds in Pmo.qp'
endif endif
A(i,j) = integral A(i,j) = integral
enddo enddo
@ -137,16 +131,16 @@ subroutine run
A = 0.d0 A = 0.d0
i = 1 i = 1
j = 1 j = 1
iunit = getunitandopen('V_mo.qp','r') iunit = getunitandopen('Vmo.qp','r')
do do
read (iunit,*,end=12) i,j, integral read (iunit,*,end=12) 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 V_mo.qp' stop 'i out of bounds in Vmo.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 V_mo.qp' stop 'j out of bounds in Vmo.qp'
endif endif
A(i,j) = integral A(i,j) = integral
enddo enddo
@ -160,7 +154,7 @@ subroutine run
! Two-electron integrals ! Two-electron integrals
iunit = getunitandopen('W_mo.qp','r') iunit = getunitandopen('Wmo.qp','r')
n_integrals=0 n_integrals=0
i = 1 i = 1
j = 1 j = 1
@ -171,25 +165,26 @@ subroutine run
read (iunit,*,end=13) i,j,k,l, integral read (iunit,*,end=13) i,j,k,l, 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 W_mo.qp' stop 'i out of bounds in Wmo.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 W_mo.qp' stop 'j out of bounds in Wmo.qp'
endif endif
if (k<0 .or. k>mo_num) then if (k<0 .or. k>mo_num) then
print *, k print *, k
stop 'k out of bounds in W_mo.qp' stop 'k out of bounds in Wmo.qp'
endif endif
if (l<0 .or. l>mo_num) then if (l<0 .or. l>mo_num) then
print *, l print *, l
stop 'l out of bounds in W_mo.qp' stop 'l out of bounds in Wmo.qp'
endif endif
n_integrals += 1 n_integrals += 1
call mo_two_e_integrals_index(i, j, k, l, buffer_i(n_integrals) ) call mo_two_e_integrals_index(i, j, k, l, buffer_i(n_integrals) )
buffer_values(n_integrals) = integral buffer_values(n_integrals) = integral
if (n_integrals == size(buffer_i)) then 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 n_integrals = 0
endif endif
enddo enddo
@ -197,15 +192,13 @@ subroutine run
close(iunit) close(iunit)
if (n_integrals > 0) then 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 endif
call map_sort(mo_integrals_map)
call map_unique(mo_integrals_map) call map_unique(mo_integrals_map)
call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints',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') call ezfio_set_mo_two_e_ints_io_mo_two_e_integrals('Read')
end end