mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-08-25 13:41:47 +02:00
working on 3idx mo ints
This commit is contained in:
parent
1c09b7dcbc
commit
727ab502c5
@ -6,69 +6,6 @@ BEGIN_PROVIDER [ integer, mo_num_per_kpt ]
|
||||
mo_num_per_kpt = mo_num/kpt_num
|
||||
END_PROVIDER
|
||||
|
||||
!BEGIN_PROVIDER [ complex*16, mo_coef_complex, (ao_num,mo_num) ]
|
||||
! implicit none
|
||||
! BEGIN_DOC
|
||||
! ! Molecular orbital coefficients on |AO| basis set
|
||||
! !
|
||||
! ! mo_coef_imag(i,j) = coefficient of the i-th |AO| on the jth |MO|
|
||||
! !
|
||||
! ! mo_label : Label characterizing the |MOs| (local, canonical, natural, etc)
|
||||
! END_DOC
|
||||
! integer :: i, j
|
||||
! double precision, allocatable :: buffer_re(:,:),buffer_im(:,:)
|
||||
! logical :: exists_re,exists_im,exists
|
||||
! PROVIDE ezfio_filename
|
||||
!
|
||||
!
|
||||
! if (mpi_master) then
|
||||
! ! Coefs
|
||||
! call ezfio_has_mo_basis_mo_coef_real(exists_re)
|
||||
! call ezfio_has_mo_basis_mo_coef_imag(exists_im)
|
||||
! exists = (exists_re.and.exists_im)
|
||||
! endif
|
||||
! IRP_IF MPI_DEBUG
|
||||
! print *, irp_here, mpi_rank
|
||||
! call MPI_BARRIER(MPI_COMM_WORLD, ierr)
|
||||
! IRP_ENDIF
|
||||
! IRP_IF MPI
|
||||
! include 'mpif.h'
|
||||
! integer :: ierr
|
||||
! call MPI_BCAST(exists, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
|
||||
! if (ierr /= MPI_SUCCESS) then
|
||||
! stop 'Unable to read mo_coef_real/imag with MPI'
|
||||
! endif
|
||||
! IRP_ENDIF
|
||||
!
|
||||
! if (exists) then
|
||||
! if (mpi_master) then
|
||||
! allocate(buffer_re(ao_num,mo_num),buffer_im(ao_num,mo_num))
|
||||
! call ezfio_get_mo_basis_mo_coef_real(buffer_re)
|
||||
! call ezfio_get_mo_basis_mo_coef_imag(buffer_im)
|
||||
! write(*,*) 'Read mo_coef_real/imag'
|
||||
! do i=1,mo_num
|
||||
! do j=1,ao_num
|
||||
! mo_coef_complex(j,i) = dcmplx(buffer_re(j,i),buffer_im(j,i))
|
||||
! enddo
|
||||
! enddo
|
||||
! deallocate(buffer_re,buffer_im)
|
||||
! endif
|
||||
! IRP_IF MPI
|
||||
! call MPI_BCAST( mo_coef_complex, mo_num*ao_num, MPI_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD, ierr)
|
||||
! if (ierr /= MPI_SUCCESS) then
|
||||
! stop 'Unable to read mo_coef_real with MPI'
|
||||
! endif
|
||||
! IRP_ENDIF
|
||||
! else
|
||||
! ! Orthonormalized AO basis
|
||||
! do i=1,mo_num
|
||||
! do j=1,ao_num
|
||||
! mo_coef_complex(j,i) = ao_ortho_canonical_coef_complex(j,i)
|
||||
! enddo
|
||||
! enddo
|
||||
! endif
|
||||
!END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ complex*16, mo_coef_complex, (ao_num,mo_num) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
@ -82,7 +19,6 @@ BEGIN_PROVIDER [ complex*16, mo_coef_complex, (ao_num,mo_num) ]
|
||||
logical :: exists
|
||||
PROVIDE ezfio_filename
|
||||
|
||||
|
||||
if (mpi_master) then
|
||||
! Coefs
|
||||
call ezfio_has_mo_basis_mo_coef_complex(exists)
|
||||
@ -121,73 +57,6 @@ BEGIN_PROVIDER [ complex*16, mo_coef_complex, (ao_num,mo_num) ]
|
||||
endif
|
||||
END_PROVIDER
|
||||
|
||||
! BEGIN_PROVIDER [ double precision, mo_coef_real, (ao_num,mo_num) ]
|
||||
!&BEGIN_PROVIDER [ double precision, mo_coef_imag, (ao_num,mo_num) ]
|
||||
!&BEGIN_PROVIDER [ complex*16, mo_coef_complex, (ao_num,mo_num) ]
|
||||
! implicit none
|
||||
! BEGIN_DOC
|
||||
! ! Molecular orbital coefficients on |AO| basis set
|
||||
! !
|
||||
! ! mo_coef_imag(i,j) = coefficient of the i-th |AO| on the jth |MO|
|
||||
! !
|
||||
! ! mo_label : Label characterizing the |MOs| (local, canonical, natural, etc)
|
||||
! END_DOC
|
||||
! integer :: i, j
|
||||
! double precision, allocatable :: buffer_re(:,:),buffer_im(:,:)
|
||||
! logical :: exists_re,exists_im,exists
|
||||
! PROVIDE ezfio_filename
|
||||
!
|
||||
!
|
||||
! if (mpi_master) then
|
||||
! ! Coefs
|
||||
! call ezfio_has_mo_basis_mo_coef_real(exists_re)
|
||||
! call ezfio_has_mo_basis_mo_coef_imag(exists_im)
|
||||
! exists = (exists_re.and.exists_im)
|
||||
! endif
|
||||
! IRP_IF MPI_DEBUG
|
||||
! print *, irp_here, mpi_rank
|
||||
! call MPI_BARRIER(MPI_COMM_WORLD, ierr)
|
||||
! IRP_ENDIF
|
||||
! IRP_IF MPI
|
||||
! include 'mpif.h'
|
||||
! integer :: ierr
|
||||
! call MPI_BCAST(exists, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
|
||||
! if (ierr /= MPI_SUCCESS) then
|
||||
! stop 'Unable to read mo_coef_real/imag with MPI'
|
||||
! endif
|
||||
! IRP_ENDIF
|
||||
!
|
||||
! if (exists) then
|
||||
! if (mpi_master) then
|
||||
! call ezfio_get_mo_basis_mo_coef_real(mo_coef_real)
|
||||
! call ezfio_get_mo_basis_mo_coef_imag(mo_coef_imag)
|
||||
! write(*,*) 'Read mo_coef_real/imag'
|
||||
! endif
|
||||
! IRP_IF MPI
|
||||
! call MPI_BCAST( mo_coef_real, mo_num*ao_num, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
|
||||
! if (ierr /= MPI_SUCCESS) then
|
||||
! stop 'Unable to read mo_coef_real with MPI'
|
||||
! endif
|
||||
! call MPI_BCAST( mo_coef_imag, mo_num*ao_num, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
|
||||
! if (ierr /= MPI_SUCCESS) then
|
||||
! stop 'Unable to read mo_coef_imag with MPI'
|
||||
! endif
|
||||
! IRP_ENDIF
|
||||
! do i=1,mo_num
|
||||
! do j=1,ao_num
|
||||
! mo_coef_complex(j,i) = dcmplx(mo_coef_real(j,i),mo_coef_imag(j,i))
|
||||
! enddo
|
||||
! enddo
|
||||
! else
|
||||
! ! Orthonormalized AO basis
|
||||
! do i=1,mo_num
|
||||
! do j=1,ao_num
|
||||
! mo_coef_complex(j,i) = ao_ortho_canonical_coef_complex(j,i)
|
||||
! enddo
|
||||
! enddo
|
||||
! endif
|
||||
!END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ complex*16, mo_coef_in_ao_ortho_basis_complex, (ao_num, mo_num) ]
|
||||
implicit none
|
||||
@ -201,6 +70,27 @@ BEGIN_PROVIDER [ complex*16, mo_coef_in_ao_ortho_basis_complex, (ao_num, mo_num)
|
||||
mo_coef_complex, size(mo_coef_complex,1), (0.d0,0.d0), &
|
||||
mo_coef_in_ao_ortho_basis_complex, size(mo_coef_in_ao_ortho_basis_complex,1))
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ complex*16, mo_coef_complex_kpts, (ao_num_per_kpt, mo_num_per_kpt, kpt_num) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! nonzero blocks of |MO| coefficients
|
||||
!
|
||||
END_DOC
|
||||
integer :: i,j,k, mo_shft, ao_shft
|
||||
mo_coef_complex_kpts = (0.d0,0.d0)
|
||||
|
||||
do k=1,kpt_num
|
||||
mo_shft = (k-1)*mo_num_per_kpt
|
||||
ao_shft = (k-1)*ao_num_per_kpt
|
||||
do i=1,mo_num_per_kpt
|
||||
do j=1,ao_num_per_kpt
|
||||
mo_coef_complex_kpts(j,i,k) = mo_coef_complex(j+ao_shft,i+mo_shft)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ complex*16, mo_coef_transp_complex, (mo_num,ao_num) ]
|
||||
|
@ -17,15 +17,9 @@ doc: Read/Write df |MO| integrals from/to disk [ Write | Read | None ]
|
||||
interface: ezfio,provider,ocaml
|
||||
default: None
|
||||
|
||||
[df_mo_integrals_real]
|
||||
[df_mo_integrals_complex]
|
||||
type: double precision
|
||||
doc: Real part of the df integrals over MOs
|
||||
size: (mo_basis.mo_num_per_kpt,mo_basis.mo_num_per_kpt,ao_two_e_ints.df_num,nuclei.kpt_pair_num)
|
||||
interface: ezfio
|
||||
|
||||
[df_mo_integrals_imag]
|
||||
type: double precision
|
||||
doc: Imaginary part of the df integrals over MOs
|
||||
size: (mo_basis.mo_num_per_kpt,mo_basis.mo_num_per_kpt,ao_two_e_ints.df_num,nuclei.kpt_pair_num)
|
||||
doc: Complex df integrals over MOs
|
||||
size: (2,mo_basis.mo_num_per_kpt,mo_basis.mo_num_per_kpt,ao_two_e_ints.df_num,nuclei.kpt_pair_num)
|
||||
interface: ezfio
|
||||
|
||||
|
@ -1,50 +1,238 @@
|
||||
BEGIN_PROVIDER [double precision, df_mo_integrals_real, (mo_num_per_kpt,mo_num_per_kpt,df_num,kpt_pair_num)]
|
||||
&BEGIN_PROVIDER [double precision, df_mo_integrals_imag, (mo_num_per_kpt,mo_num_per_kpt,df_num,kpt_pair_num)]
|
||||
&BEGIN_PROVIDER [complex*16, df_mo_integrals_complex, (mo_num_per_kpt,mo_num_per_kpt,df_num,kpt_pair_num)]
|
||||
BEGIN_PROVIDER [complex*16, df_mo_integrals_complex, (mo_num_per_kpt,mo_num_per_kpt,df_num,kpt_pair_num)]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! df AO integrals
|
||||
! df MO integrals
|
||||
END_DOC
|
||||
integer :: i,j,k,l
|
||||
|
||||
if (read_df_mo_integrals) then
|
||||
df_mo_integrals_real = 0.d0
|
||||
df_mo_integrals_imag = 0.d0
|
||||
call ezfio_get_mo_two_e_ints_df_mo_integrals_real(df_mo_integrals_real)
|
||||
call ezfio_get_mo_two_e_ints_df_mo_integrals_imag(df_mo_integrals_imag)
|
||||
print *, 'df AO integrals read from disk'
|
||||
do l=1,kpt_pair_num
|
||||
do k=1,df_num
|
||||
do j=1,mo_num_per_kpt
|
||||
do i=1,mo_num_per_kpt
|
||||
df_mo_integrals_complex(i,j,k,l) = dcmplx(df_mo_integrals_real(i,j,k,l), &
|
||||
df_mo_integrals_imag(i,j,k,l))
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
call ezfio_get_mo_two_e_ints_df_mo_integrals_complex(df_mo_integrals_complex)
|
||||
print *, 'df MO integrals read from disk'
|
||||
else
|
||||
call df_mo_from_df_ao(df_mo_integrals_complex,df_ao_integrals_complex,mo_num_per_kpt,ao_num_per_kpt,df_num,kpt_pair_num)
|
||||
endif
|
||||
|
||||
if (write_df_mo_integrals) then
|
||||
do l=1,kpt_pair_num
|
||||
do k=1,df_num
|
||||
do j=1,mo_num_per_kpt
|
||||
do i=1,mo_num_per_kpt
|
||||
df_mo_integrals_real(i,j,k,l) = dble(df_mo_integrals_complex(i,j,k,l))
|
||||
df_mo_integrals_imag(i,j,k,l) = dimag(df_mo_integrals_complex(i,j,k,l))
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
call ezfio_set_mo_two_e_ints_df_mo_integrals_real(df_mo_integrals_real)
|
||||
call ezfio_set_mo_two_e_ints_df_mo_integrals_imag(df_mo_integrals_imag)
|
||||
print *, 'df AO integrals written to disk'
|
||||
call ezfio_set_mo_two_e_ints_df_mo_integrals_complex(df_mo_integrals_complex)
|
||||
print *, 'df MO integrals written to disk'
|
||||
endif
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
subroutine mo_map_fill_from_df
|
||||
use map_module
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! fill mo bielec integral map using 3-index df integrals
|
||||
END_DOC
|
||||
|
||||
integer :: i,k,j,l
|
||||
integer :: ki,kk,kj,kl
|
||||
integer :: ii,ik,ij,il
|
||||
integer :: kikk2,kjkl2,jl2,ik2
|
||||
integer :: i_mo,j_mo,i_df
|
||||
|
||||
complex*16,allocatable :: ints_ik(:,:,:), ints_jl(:,:,:), ints_ikjl(:,:,:,:)
|
||||
|
||||
complex*16 :: integral
|
||||
integer :: n_integrals_1, n_integrals_2
|
||||
integer :: size_buffer
|
||||
integer(key_kind),allocatable :: buffer_i_1(:), buffer_i_2(:)
|
||||
real(integral_kind),allocatable :: buffer_values_1(:), buffer_values_2(:)
|
||||
double precision :: tmp_re,tmp_im
|
||||
integer :: mo_num_kpt_2
|
||||
|
||||
double precision :: cpu_1, cpu_2, wall_1, wall_2, wall_0
|
||||
double precision :: map_mb
|
||||
|
||||
logical :: use_map1
|
||||
integer(keY_kind) :: idx_tmp
|
||||
double precision :: sign
|
||||
|
||||
mo_num_kpt_2 = mo_num_per_kpt * mo_num_per_kpt
|
||||
|
||||
size_buffer = min(mo_num_per_kpt*mo_num_per_kpt*mo_num_per_kpt,16000000)
|
||||
print*, 'Providing the mo_bielec integrals from 3-index df integrals'
|
||||
call write_time(6)
|
||||
! call ezfio_set_integrals_bielec_disk_access_mo_integrals('Write')
|
||||
! TOUCH read_mo_integrals read_ao_integrals write_mo_integrals write_ao_integrals
|
||||
|
||||
call wall_time(wall_1)
|
||||
call cpu_time(cpu_1)
|
||||
|
||||
allocate( ints_jl(mo_num_per_kpt,mo_num_per_kpt,df_num))
|
||||
|
||||
wall_0 = wall_1
|
||||
do kl=1, kpt_num
|
||||
do kj=1, kl
|
||||
call idx2_tri_int(kj,kl,kjkl2)
|
||||
if (kj < kl) then
|
||||
do i_mo=1,mo_num_per_kpt
|
||||
do j_mo=1,mo_num_per_kpt
|
||||
do i_df=1,df_num
|
||||
ints_jl(i_mo,j_mo,i_df) = dconjg(df_mo_integrals_complex(j_mo,i_mo,i_df,kjkl2))
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
else
|
||||
ints_jl = df_mo_integrals_complex(:,:,:,kjkl2)
|
||||
endif
|
||||
|
||||
!$OMP PARALLEL PRIVATE(i,k,j,l,ki,kk,ii,ik,ij,il,kikk2,jl2,ik2, &
|
||||
!$OMP ints_ik, ints_ikjl, i_mo, j_mo, i_df, &
|
||||
!$OMP n_integrals_1, buffer_i_1, buffer_values_1, &
|
||||
!$OMP n_integrals_2, buffer_i_2, buffer_values_2, &
|
||||
!$OMP idx_tmp, tmp_re, tmp_im, integral,sign,use_map1) &
|
||||
!$OMP DEFAULT(NONE) &
|
||||
!$OMP SHARED(size_buffer, kpt_num, df_num, mo_num_per_kpt, mo_num_kpt_2, &
|
||||
!$OMP kl,kj,kjkl2,ints_jl, &
|
||||
!$OMP kconserv, df_mo_integrals_complex, mo_integrals_threshold, mo_integrals_map, mo_integrals_map_2)
|
||||
|
||||
allocate( &
|
||||
ints_ik(mo_num_per_kpt,mo_num_per_kpt,df_num), &
|
||||
ints_ikjl(mo_num_per_kpt,mo_num_per_kpt,mo_num_per_kpt,mo_num_per_kpt), &
|
||||
buffer_i_1(size_buffer), &
|
||||
buffer_i_2(size_buffer), &
|
||||
buffer_values_1(size_buffer), &
|
||||
buffer_values_2(size_buffer) &
|
||||
)
|
||||
|
||||
!$OMP DO SCHEDULE(guided)
|
||||
do kk=1,kl
|
||||
ki=kconserv(kl,kk,kj)
|
||||
if (ki>kl) cycle
|
||||
! if ((kl == kj) .and. (ki > kk)) cycle
|
||||
call idx2_tri_int(ki,kk,kikk2)
|
||||
! if (kikk2 > kjkl2) cycle
|
||||
if (ki < kk) then
|
||||
do i_mo=1,mo_num_per_kpt
|
||||
do j_mo=1,mo_num_per_kpt
|
||||
do i_df=1,df_num
|
||||
ints_ik(i_mo,j_mo,i_df) = dconjg(df_mo_integrals_complex(j_mo,i_mo,i_df,kikk2))
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
! ints_ik = conjg(reshape(df_mo_integral_array(:,:,:,kikk2),(/mo_num_per_kpt,mo_num_per_kpt,df_num/),order=(/2,1,3/)))
|
||||
else
|
||||
ints_ik = df_mo_integrals_complex(:,:,:,kikk2)
|
||||
endif
|
||||
|
||||
call zgemm('N','T', mo_num_kpt_2, mo_num_kpt_2, df_num, &
|
||||
(1.d0,0.d0), ints_ik, mo_num_kpt_2, &
|
||||
ints_jl, mo_num_kpt_2, &
|
||||
(0.d0,0.d0), ints_ikjl, mo_num_kpt_2)
|
||||
|
||||
n_integrals_1=0
|
||||
n_integrals_2=0
|
||||
do il=1,mo_num_per_kpt
|
||||
l=il+(kl-1)*mo_num_per_kpt
|
||||
do ij=1,mo_num_per_kpt
|
||||
j=ij+(kj-1)*mo_num_per_kpt
|
||||
if (j>l) exit
|
||||
call idx2_tri_int(j,l,jl2)
|
||||
do ik=1,mo_num_per_kpt
|
||||
k=ik+(kk-1)*mo_num_per_kpt
|
||||
if (k>l) exit
|
||||
do ii=1,mo_num_per_kpt
|
||||
i=ii+(ki-1)*mo_num_per_kpt
|
||||
if ((j==l) .and. (i>k)) exit
|
||||
call idx2_tri_int(i,k,ik2)
|
||||
if (ik2 > jl2) exit
|
||||
integral = ints_ikjl(ii,ik,ij,il)
|
||||
! print*,i,k,j,l,real(integral),imag(integral)
|
||||
if (cdabs(integral) < mo_integrals_threshold) then
|
||||
cycle
|
||||
endif
|
||||
call ao_two_e_integral_complex_map_idx_sign(i,j,k,l,use_map1,idx_tmp,sign)
|
||||
tmp_re = dble(integral)
|
||||
tmp_im = dimag(integral)
|
||||
if (use_map1) then
|
||||
n_integrals_1 += 1
|
||||
buffer_i_1(n_integrals_1)=idx_tmp
|
||||
buffer_values_1(n_integrals_1)=tmp_re
|
||||
if (sign.ne.0.d0) then
|
||||
n_integrals_1 += 1
|
||||
buffer_i_1(n_integrals_1)=idx_tmp+1
|
||||
buffer_values_1(n_integrals_1)=tmp_im*sign
|
||||
endif
|
||||
if (n_integrals_1 >= size(buffer_i_1)-1) then
|
||||
call map_append(mo_integrals_map, buffer_i_1, buffer_values_1, n_integrals_1)
|
||||
!call insert_into_ao_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1)
|
||||
n_integrals_1 = 0
|
||||
endif
|
||||
else
|
||||
n_integrals_2 += 1
|
||||
buffer_i_2(n_integrals_2)=idx_tmp
|
||||
buffer_values_2(n_integrals_2)=tmp_re
|
||||
if (sign.ne.0.d0) then
|
||||
n_integrals_2 += 1
|
||||
buffer_i_2(n_integrals_2)=idx_tmp+1
|
||||
buffer_values_2(n_integrals_2)=tmp_im*sign
|
||||
endif
|
||||
if (n_integrals_2 >= size(buffer_i_2)-1) then
|
||||
call map_append(mo_integrals_map_2, buffer_i_2, buffer_values_2, n_integrals_2)
|
||||
!call insert_into_ao_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_2)
|
||||
n_integrals_2 = 0
|
||||
endif
|
||||
endif
|
||||
|
||||
enddo !ii
|
||||
enddo !ik
|
||||
enddo !ij
|
||||
enddo !il
|
||||
|
||||
if (n_integrals_1 > 0) then
|
||||
call map_append(mo_integrals_map, buffer_i_1, buffer_values_1, n_integrals_1)
|
||||
!call insert_into_ao_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1)
|
||||
endif
|
||||
if (n_integrals_2 > 0) then
|
||||
call map_append(mo_integrals_map_2, buffer_i_2, buffer_values_2, n_integrals_2)
|
||||
!call insert_into_ao_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_2)
|
||||
endif
|
||||
enddo !kk
|
||||
!$OMP END DO NOWAIT
|
||||
deallocate( &
|
||||
ints_ik, &
|
||||
ints_ikjl, &
|
||||
buffer_i_1, &
|
||||
buffer_i_2, &
|
||||
buffer_values_1, &
|
||||
buffer_values_2 &
|
||||
)
|
||||
!$OMP END PARALLEL
|
||||
enddo !kj
|
||||
call wall_time(wall_2)
|
||||
if (wall_2 - wall_0 > 1.d0) then
|
||||
wall_0 = wall_2
|
||||
print*, 100.*float(kl)/float(kpt_num), '% in ', &
|
||||
wall_2-wall_1,'s',map_mb(mo_integrals_map),'+',map_mb(mo_integrals_map_2),'MB'
|
||||
endif
|
||||
|
||||
enddo !kl
|
||||
deallocate( ints_jl )
|
||||
|
||||
call map_sort(mo_integrals_map)
|
||||
call map_unique(mo_integrals_map)
|
||||
call map_sort(mo_integrals_map_2)
|
||||
call map_unique(mo_integrals_map_2)
|
||||
!call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints_complex_1',mo_integrals_map)
|
||||
!call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints_complex_2',mo_integrals_map_2)
|
||||
!call ezfio_set_mo_two_e_ints_io_mo_two_e_integrals('Read')
|
||||
|
||||
call wall_time(wall_2)
|
||||
call cpu_time(cpu_2)
|
||||
|
||||
integer*8 :: get_mo_map_size, mo_map_size
|
||||
mo_map_size = get_mo_map_size()
|
||||
|
||||
print*,'MO integrals provided:'
|
||||
print*,' Size of MO map ', map_mb(mo_integrals_map),'+',map_mb(mo_integrals_map_2),'MB'
|
||||
print*,' Number of MO integrals: ', mo_map_size
|
||||
print*,' cpu time :',cpu_2 - cpu_1, 's'
|
||||
print*,' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1), ')'
|
||||
|
||||
end subroutine mo_map_fill_from_df
|
||||
|
||||
subroutine df_mo_from_df_ao(df_mo,df_ao,n_mo,n_ao,n_df,n_k_pairs)
|
||||
use map_module
|
||||
implicit none
|
||||
@ -70,9 +258,9 @@ subroutine df_mo_from_df_ao(df_mo,df_ao,n_mo,n_ao,n_df,n_k_pairs)
|
||||
)
|
||||
|
||||
do kl=1, kpt_num
|
||||
coef_l = mo_coef_kpts(:,:,kl)
|
||||
coef_l = mo_coef_complex_kpts(:,:,kl)
|
||||
do kj=1, kl
|
||||
coef_j = mo_coef_kpts(:,:,kj)
|
||||
coef_j = mo_coef_complex_kpts(:,:,kj)
|
||||
kjkl2 = kj+shiftr(kl*kl-kl,1)
|
||||
do mu=1, df_num
|
||||
ints_jl = df_ao(:,:,mu,kjkl2)
|
||||
|
@ -41,6 +41,10 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ]
|
||||
call map_load_from_disk(trim(ezfio_filename)//'/work/mo_ints_complex_2',mo_integrals_map_2)
|
||||
print*, 'MO integrals provided (periodic)'
|
||||
return
|
||||
else if (read_df_mo_integrals) then
|
||||
PROVIDE df_mo_integrals_complex
|
||||
call mo_map_fill_from_df
|
||||
return
|
||||
else
|
||||
PROVIDE ao_two_e_integrals_in_map
|
||||
endif
|
||||
|
@ -36,6 +36,9 @@ def idx4(i,j,k,l):
|
||||
def stri4z(i,j,k,l,zr,zi):
|
||||
return (4*'{:5d}'+2*'{:25.16e}').format(i,j,k,l,zr,zi)
|
||||
|
||||
def stri2z(i,j,zr,zi):
|
||||
return (2*'{:5d}'+2*'{:25.16e}').format(i,j,zr,zi)
|
||||
|
||||
def strijklikjli4z(i,j,k,l,zr,zi):
|
||||
return ('{:10d}'+ 2*'{:8d}'+4*'{:5d}'+2*'{:25.16e}').format(idx4(i,j,k,l),idx2_tri((i-1,k-1))+1,idx2_tri((j-1,l-1))+1,i,j,k,l,zr,zi)
|
||||
|
||||
@ -322,7 +325,7 @@ def pyscf2QP(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8,
|
||||
for j in range(i,nao):
|
||||
int_ij = intval_kpts_ao[ik,i,j]
|
||||
if abs(int_ij) > thresh:
|
||||
outfile.write('%s %s %s %s\n' % (i+shift, j+shift, int_ij.real, int_ij.imag))
|
||||
outfile.write(stri2z(i+shift, j+shift, int_ij.real, int_ij.imag)+'\n')
|
||||
if print_mo_ints_mono:
|
||||
intval_kpts_mo = np.einsum('kim,kij,kjn->kmn',mo_k.conj(),intval_kpts_ao,mo_k)
|
||||
with open('%s_mo_complex' % name,'w') as outfile:
|
||||
@ -332,7 +335,7 @@ def pyscf2QP(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8,
|
||||
for j in range(i,nmo):
|
||||
int_ij = intval_kpts_mo[ik,i,j]
|
||||
if abs(int_ij) > thresh:
|
||||
outfile.write('%s %s %s %s\n' % (i+shift, j+shift, int_ij.real, int_ij.imag))
|
||||
outfile.write(stri2z(i+shift, j+shift, int_ij.real, int_ij.imag)+'\n')
|
||||
|
||||
|
||||
# ___ _
|
||||
@ -383,7 +386,7 @@ def pyscf2QP(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8,
|
||||
for i,i0 in enumerate(dfbasfunc):
|
||||
for j,v in enumerate(i0):
|
||||
if (abs(v) > bielec_int_threshold):
|
||||
outfile.write('%s %s %s %s %s %s\n' % (i+1,j+1,iaux+1,k+1,v.real,v.imag))
|
||||
outfile.write(stri4z(i+1,j+1,iaux+1,k+1,v.real,v.imag)+'\n')
|
||||
|
||||
if print_mo_ints_df:
|
||||
kpair_list=[]
|
||||
@ -400,28 +403,9 @@ def pyscf2QP(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8,
|
||||
for i,i0 in enumerate(dfbasfunc):
|
||||
for j,v in enumerate(i0):
|
||||
if (abs(v) > bielec_int_threshold):
|
||||
outfile.write('%s %s %s %s %s %s\n' % (i+1,j+1,iaux+1,k+1,v.real,v.imag))
|
||||
outfile.write(stri4z(i+1,j+1,iaux+1,k+1,v.real,v.imag)+'\n')
|
||||
|
||||
|
||||
|
||||
|
||||
# eri_4d_ao = np.zeros((Nk,nao,Nk,nao,Nk,nao,Nk,nao), dtype=np.complex)
|
||||
# for d, kd in enumerate(kpts):
|
||||
# for c, kc in enumerate(kpts):
|
||||
# if c > d: break
|
||||
# idx2_cd = idx2_tri(c,d)
|
||||
# for b, kb in enumerate(kpts):
|
||||
# if b > d: break
|
||||
# a = kconserv[b,c,d]
|
||||
# if idx2_tri(a,b) > idx2_cd: continue
|
||||
# if ((c==d) and (a>b)): continue
|
||||
# ka = kpts[a]
|
||||
# v = mf.with_df.get_ao_eri(kpts=[ka,kb,kc,kd],compact=False).reshape((nao,)*4)
|
||||
# v *= 1./Nk
|
||||
# eri_4d_ao[a,:,b,:,c,:,d] = v
|
||||
#
|
||||
# eri_4d_ao = eri_4d_ao.reshape([Nk*nao]*4)
|
||||
|
||||
if (print_ao_ints_bi):
|
||||
print_ao_bi(mf,kconserv,'bielec_ao_complex',bielec_int_threshold)
|
||||
if (print_mo_ints_bi):
|
||||
@ -446,7 +430,6 @@ def print_mo_bi(mf,kconserv=None,outfilename='W.mo.qp',cas_idx=None,bielec_int_t
|
||||
kconserv = tools.get_kconserv(cell, kpts)
|
||||
|
||||
with open(outfilename,'w') as outfile:
|
||||
pass
|
||||
for d, kd in enumerate(kpts):
|
||||
for c, kc in enumerate(kpts):
|
||||
if c > d: break
|
||||
@ -458,7 +441,6 @@ def print_mo_bi(mf,kconserv=None,outfilename='W.mo.qp',cas_idx=None,bielec_int_t
|
||||
#if idx2_tri((a,b)) > idx2_cd: continue
|
||||
#if ((c==d) and (a>b)): continue
|
||||
ka = kpts[a]
|
||||
with open(outfilename,'a') as outfile:
|
||||
eri_4d_mo_kpt = mf.with_df.ao2mo([mo_k[a], mo_k[b], mo_k[c], mo_k[d]],
|
||||
[ka,kb,kc,kd],compact=False).reshape((nmo,)*4)
|
||||
eri_4d_mo_kpt *= 1./Nk
|
||||
@ -477,7 +459,8 @@ def print_mo_bi(mf,kconserv=None,outfilename='W.mo.qp',cas_idx=None,bielec_int_t
|
||||
if ((jj==ll) and (ii>kk)): break
|
||||
v=eri_4d_mo_kpt[i,k,j,l]
|
||||
if (abs(v) > bielec_int_threshold):
|
||||
outfile.write(stri4z(ii+1,jj+1,kk+1,ll+1,v.real,v.imag)+'\n')
|
||||
outfile.write(stri4z(ii+1,jj+1,kk+1,ll+1,
|
||||
v.real,v.imag)+'\n')
|
||||
|
||||
|
||||
def print_ao_bi(mf,kconserv=None,outfilename='W.ao.qp',bielec_int_threshold = 1E-8):
|
||||
@ -492,7 +475,6 @@ def print_ao_bi(mf,kconserv=None,outfilename='W.ao.qp',bielec_int_threshold = 1E
|
||||
kconserv = tools.get_kconserv(cell, kpts)
|
||||
|
||||
with open(outfilename,'w') as outfile:
|
||||
pass
|
||||
for d, kd in enumerate(kpts):
|
||||
for c, kc in enumerate(kpts):
|
||||
if c > d: break
|
||||
@ -505,8 +487,8 @@ def print_ao_bi(mf,kconserv=None,outfilename='W.ao.qp',bielec_int_threshold = 1E
|
||||
#if ((c==d) and (a>b)): continue
|
||||
ka = kpts[a]
|
||||
|
||||
with open(outfilename,'a') as outfile:
|
||||
eri_4d_ao_kpt = mf.with_df.get_ao_eri(kpts=[ka,kb,kc,kd],compact=False).reshape((nao,)*4)
|
||||
eri_4d_ao_kpt = mf.with_df.get_ao_eri(kpts=[ka,kb,kc,kd],
|
||||
compact=False).reshape((nao,)*4)
|
||||
eri_4d_ao_kpt *= 1./Nk
|
||||
for l in range(nao):
|
||||
ll=l+d*nao
|
||||
@ -523,7 +505,8 @@ def print_ao_bi(mf,kconserv=None,outfilename='W.ao.qp',bielec_int_threshold = 1E
|
||||
if ((jj==ll) and (ii>kk)): break
|
||||
v=eri_4d_ao_kpt[i,k,j,l]
|
||||
if (abs(v) > bielec_int_threshold):
|
||||
outfile.write(stri4z(ii+1,jj+1,kk+1,ll+1,v.real,v.imag)+'\n')
|
||||
outfile.write(stri4z(ii+1,jj+1,kk+1,ll+1,
|
||||
v.real,v.imag)+'\n')
|
||||
|
||||
|
||||
|
||||
@ -640,7 +623,7 @@ def pyscf2QP2(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8,
|
||||
for j in range(nmo):
|
||||
cij = c_kpts[ik,i,j]
|
||||
if abs(cij) > mo_coef_threshold:
|
||||
outfile.write('%s %s %s %s\n' % (i+shift1, j+shift2, cij.real, cij.imag))
|
||||
outfile.write(stri2z(i+shift1, j+shift2, cij.real, cij.imag)+'\n')
|
||||
|
||||
# ___
|
||||
# | ._ _|_ _ _ ._ _. | _ |\/| _ ._ _
|
||||
@ -679,7 +662,7 @@ def pyscf2QP2(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8,
|
||||
for j in range(i,nao):
|
||||
int_ij = intval_kpts_ao[ik,i,j]
|
||||
if abs(int_ij) > thresh:
|
||||
outfile.write('%s %s %s %s\n' % (i+shift, j+shift, int_ij.real, int_ij.imag))
|
||||
outfile.write(stri2z(i+shift, j+shift, int_ij.real, int_ij.imag)+'\n')
|
||||
if print_mo_ints_mono:
|
||||
intval_kpts_mo = np.einsum('kim,kij,kjn->kmn',mo_k.conj(),intval_kpts_ao,mo_k)
|
||||
with open('%s_mo.qp' % name,'w') as outfile:
|
||||
@ -689,7 +672,7 @@ def pyscf2QP2(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8,
|
||||
for j in range(i,nmo):
|
||||
int_ij = intval_kpts_mo[ik,i,j]
|
||||
if abs(int_ij) > thresh:
|
||||
outfile.write('%s %s %s %s\n' % (i+shift, j+shift, int_ij.real, int_ij.imag))
|
||||
outfile.write(stri2z(i+shift, j+shift, int_ij.real, int_ij.imag)+'\n')
|
||||
|
||||
|
||||
# ___ _
|
||||
@ -749,7 +732,7 @@ def pyscf2QP2(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8,
|
||||
for i,i0 in enumerate(dfbasfunc):
|
||||
for j,v in enumerate(i0):
|
||||
if (abs(v) > bielec_int_threshold):
|
||||
outfile.write('%s %s %s %s %s %s\n' % (i+1,j+1,iaux+1,k+1,v.real,v.imag))
|
||||
outfile.write(stri4z(i+1,j+1,iaux+1,k+1,v.real,v.imag)+'\n')
|
||||
df_ao_tmp[i,j,iaux,k]=v
|
||||
|
||||
qph5.create_dataset('ao_two_e_ints/df_ao_integrals_real',data=df_ao_tmp.real)
|
||||
@ -771,92 +754,15 @@ def pyscf2QP2(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8,
|
||||
for i,i0 in enumerate(dfbasfunc):
|
||||
for j,v in enumerate(i0):
|
||||
if (abs(v) > bielec_int_threshold):
|
||||
outfile.write('%s %s %s %s %s %s\n' % (i+1,j+1,iaux+1,k+1,v.real,v.imag))
|
||||
outfile.write(stri4z(i+1,j+1,iaux+1,k+1,v.real,v.imag)+'\n')
|
||||
df_mo_tmp[i,j,iaux,k]=v
|
||||
qph5.create_dataset('mo_two_e_ints/df_mo_integrals_real',data=df_mo_tmp.real)
|
||||
qph5.create_dataset('mo_two_e_ints/df_mo_integrals_imag',data=df_mo_tmp.imag)
|
||||
|
||||
|
||||
|
||||
# eri_4d_ao = np.zeros((Nk,nao,Nk,nao,Nk,nao,Nk,nao), dtype=np.complex)
|
||||
# for d, kd in enumerate(kpts):
|
||||
# for c, kc in enumerate(kpts):
|
||||
# if c > d: break
|
||||
# idx2_cd = idx2_tri(c,d)
|
||||
# for b, kb in enumerate(kpts):
|
||||
# if b > d: break
|
||||
# a = kconserv[b,c,d]
|
||||
# if idx2_tri(a,b) > idx2_cd: continue
|
||||
# if ((c==d) and (a>b)): continue
|
||||
# ka = kpts[a]
|
||||
# v = mf.with_df.get_ao_eri(kpts=[ka,kb,kc,kd],compact=False).reshape((nao,)*4)
|
||||
# v *= 1./Nk
|
||||
# eri_4d_ao[a,:,b,:,c,:,d] = v
|
||||
#
|
||||
# eri_4d_ao = eri_4d_ao.reshape([Nk*nao]*4)
|
||||
|
||||
|
||||
if (print_ao_ints_bi or print_mo_ints_bi):
|
||||
if print_ao_ints_bi:
|
||||
with open('W.qp','w') as outfile:
|
||||
pass
|
||||
if print_mo_ints_bi:
|
||||
with open('W_mo.qp','w') as outfile:
|
||||
pass
|
||||
for d, kd in enumerate(kpts):
|
||||
for c, kc in enumerate(kpts):
|
||||
if c > d: break
|
||||
idx2_cd = idx2_tri((c,d))
|
||||
for b, kb in enumerate(kpts):
|
||||
if b > d: break
|
||||
a = kconserv[b,c,d]
|
||||
#if idx2_tri((a,b)) > idx2_cd: continue
|
||||
if a > d: continue
|
||||
#if ((c==d) and (a>b)): continue
|
||||
ka = kpts[a]
|
||||
|
||||
if print_ao_ints_bi:
|
||||
with open('W.qp','a') as outfile:
|
||||
eri_4d_ao_kpt = mf.with_df.get_ao_eri(kpts=[ka,kb,kc,kd],compact=False).reshape((nao,)*4)
|
||||
eri_4d_ao_kpt *= 1./Nk
|
||||
for l in range(nao):
|
||||
ll=l+d*nao
|
||||
for j in range(nao):
|
||||
jj=j+c*nao
|
||||
if jj>ll: break
|
||||
idx2_jjll = idx2_tri((jj,ll))
|
||||
for k in range(nao):
|
||||
kk=k+b*nao
|
||||
if kk>ll: break
|
||||
for i in range(nao):
|
||||
ii=i+a*nao
|
||||
if idx2_tri((ii,kk)) > idx2_jjll: break
|
||||
if ((jj==ll) and (ii>kk)): break
|
||||
v=eri_4d_ao_kpt[i,k,j,l]
|
||||
if (abs(v) > bielec_int_threshold):
|
||||
outfile.write('%s %s %s %s %s %s\n' % (ii+1,jj+1,kk+1,ll+1,v.real,v.imag))
|
||||
|
||||
if print_mo_ints_bi:
|
||||
with open('W_mo.qp','a') as outfile:
|
||||
eri_4d_mo_kpt = mf.with_df.ao2mo([mo_k[a], mo_k[b], mo_k[c], mo_k[d]],
|
||||
[ka,kb,kc,kd],compact=False).reshape((nmo,)*4)
|
||||
eri_4d_mo_kpt *= 1./Nk
|
||||
for l in range(nmo):
|
||||
ll=l+d*nmo
|
||||
for j in range(nmo):
|
||||
jj=j+c*nmo
|
||||
if jj>ll: break
|
||||
idx2_jjll = idx2_tri((jj,ll))
|
||||
for k in range(nmo):
|
||||
kk=k+b*nmo
|
||||
if kk>ll: break
|
||||
for i in range(nmo):
|
||||
ii=i+a*nmo
|
||||
if idx2_tri((ii,kk)) > idx2_jjll: break
|
||||
if ((jj==ll) and (ii>kk)): break
|
||||
v=eri_4d_mo_kpt[i,k,j,l]
|
||||
if (abs(v) > bielec_int_threshold):
|
||||
outfile.write('%s %s %s %s %s %s\n' % (ii+1,jj+1,kk+1,ll+1,v.real,v.imag))
|
||||
if (print_ao_ints_bi):
|
||||
print_ao_bi(mf,kconserv,'W.qp',bielec_int_threshold)
|
||||
if (print_mo_ints_bi):
|
||||
print_mo_bi(mf,kconserv,'W.mo.qp',cas_idx,bielec_int_threshold)
|
||||
|
||||
|
||||
|
||||
|
@ -1,2 +1,3 @@
|
||||
ao_two_e_ints
|
||||
ao_one_e_ints
|
||||
mo_two_e_ints
|
||||
|
@ -21,9 +21,6 @@ ezfio.set_nuclei_kpt_num(kpt_num)
|
||||
kpt_pair_num = (kpt_num*kpt_num + kpt_num)//2
|
||||
ezfio.set_nuclei_kpt_pair_num(kpt_pair_num)
|
||||
|
||||
# should this be in ao_basis? ao_two_e_ints?
|
||||
df_num = qph5['ao_two_e_ints'].attrs['df_num']
|
||||
ezfio.set_ao_two_e_ints_df_num(df_num)
|
||||
|
||||
# these are totals (kpt_num * num_per_kpt)
|
||||
# need to change if we want to truncate orbital space within pyscf
|
||||
@ -139,14 +136,25 @@ ezfio.set_ao_one_e_ints_io_ao_integrals_kinetic('Read')
|
||||
ezfio.set_ao_one_e_ints_io_ao_integrals_overlap('Read')
|
||||
ezfio.set_ao_one_e_ints_io_ao_integrals_n_e('Read')
|
||||
|
||||
# should this be in ao_basis? ao_two_e_ints?
|
||||
if 'ao_two_e_ints' in qph5.keys():
|
||||
df_num = qph5['ao_two_e_ints'].attrs['df_num']
|
||||
ezfio.set_ao_two_e_ints_df_num(df_num)
|
||||
if 'df_ao_integrals_real' in qph5['ao_two_e_ints'].keys():
|
||||
dfao_re0=qph5['ao_two_e_ints/df_ao_integrals_real'][()].transpose((3,2,1,0))
|
||||
dfao_im0=qph5['ao_two_e_ints/df_ao_integrals_imag'][()].transpose((3,2,1,0))
|
||||
#ezfio.set_ao_two_e_ints_df_ao_integrals_real(dfao_re.tolist())
|
||||
#ezfio.set_ao_two_e_ints_df_ao_integrals_imag(dfao_im.tolist())
|
||||
dfao_cmplx0 = np.stack((dfao_re0,dfao_im0),axis=-1).tolist()
|
||||
ezfio.set_ao_two_e_ints_df_ao_integrals_complex(dfao_cmplx0)
|
||||
ezfio.set_ao_two_e_ints_io_df_ao_integrals('Read')
|
||||
|
||||
if 'mo_two_e_ints' in qph5.keys():
|
||||
df_num = qph5['ao_two_e_ints'].attrs['df_num']
|
||||
dfmo_re0=qph5['mo_two_e_ints/df_mo_integrals_real'][()].transpose((3,2,1,0))
|
||||
dfmo_im0=qph5['mo_two_e_ints/df_mo_integrals_imag'][()].transpose((3,2,1,0))
|
||||
dfmo_cmplx0 = np.stack((dfmo_re0,dfmo_im0),axis=-1).tolist()
|
||||
ezfio.set_mo_two_e_ints_df_mo_integrals_complex(dfmo_cmplx0)
|
||||
ezfio.set_mo_two_e_ints_io_df_mo_integrals('Read')
|
||||
|
||||
|
||||
#TODO: add check and only do this if ints exist
|
||||
#dfmo_re=qph5['mo_two_e_ints/df_mo_integrals_real'][()].transpose((3,2,1,0)).tolist()
|
||||
|
Loading…
Reference in New Issue
Block a user