mirror of
https://gitlab.com/scemama/qp_plugins_scemama.git
synced 2024-06-02 03:15:25 +02:00
94 lines
2.1 KiB
Fortran
94 lines
2.1 KiB
Fortran
program export_integrals_ao
|
|
call run
|
|
end
|
|
|
|
subroutine write_2d(name,A)
|
|
implicit none
|
|
character*(*), intent(in) :: name
|
|
double precision, intent(in) :: A(ao_num,ao_num)
|
|
integer :: i, j
|
|
integer :: iunit
|
|
integer :: getunitandopen
|
|
|
|
iunit = getunitandopen(name,'w')
|
|
do j=1,ao_num
|
|
do i=1,ao_num
|
|
if (A(i,j) /= 0.d0) then
|
|
write (iunit,*) i,j, A(i,j)
|
|
endif
|
|
enddo
|
|
enddo
|
|
close(iunit)
|
|
end
|
|
|
|
subroutine run
|
|
use map_module
|
|
implicit none
|
|
BEGIN_DOC
|
|
! Program to import integrals in the AO basis.
|
|
!
|
|
! one-electron integrals : format is : i j value
|
|
! two-electron integrals : format is : i j k l value
|
|
! Dirac's notation is used : <ij|kl> is <r1 r2|r1 r2>
|
|
!
|
|
! These files are read:
|
|
!
|
|
! E.qp : Contains the nuclear repulsion energy
|
|
!
|
|
! T.qp : kinetic energy integrals
|
|
!
|
|
! S.qp : overlap matrix
|
|
!
|
|
! P.qp : pseudopotential integrals
|
|
!
|
|
! V.qp : electron-nucleus potential
|
|
!
|
|
! W.qp : electron repulsion integrals
|
|
!
|
|
END_DOC
|
|
|
|
integer :: iunit
|
|
integer :: getunitandopen
|
|
|
|
integer :: i,j,k,l
|
|
double precision :: integral
|
|
double precision, allocatable :: A(:,:)
|
|
|
|
integer :: n_integrals
|
|
integer(key_kind) :: idx
|
|
|
|
double precision, external :: get_ao_two_e_integral
|
|
|
|
allocate (A(ao_num,ao_num))
|
|
call ezfio_set_nuclei_nuclear_repulsion(nuclear_repulsion)
|
|
|
|
call write_2d('T.qp',ao_kinetic_integrals)
|
|
call write_2d('S.qp',ao_overlap)
|
|
call write_2d('P.qp',ao_pseudo_integrals)
|
|
call write_2d('V.qp',ao_integrals_n_e)
|
|
|
|
iunit = getunitandopen('E.qp','w')
|
|
write(iunit,*) nuclear_repulsion
|
|
close(iunit)
|
|
|
|
iunit = getunitandopen('W.qp','w')
|
|
|
|
do l=1,ao_num
|
|
do k=1,ao_num
|
|
do j=l,ao_num
|
|
do i=k,ao_num
|
|
if (i==j .and. k<l) cycle
|
|
if (i>=j) then
|
|
integral = get_ao_two_e_integral(i,j,k,l,ao_integrals_map)
|
|
if (integral /= 0.d0) then
|
|
write (iunit,'(4(I5,2X),E22.15)') i,j,k,l, integral
|
|
endif
|
|
endif
|
|
enddo
|
|
enddo
|
|
enddo
|
|
enddo
|
|
close(iunit)
|
|
|
|
end
|