From 727ab502c55d49067a5b66795ca47f37042fb900 Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Tue, 18 Feb 2020 18:32:47 -0600 Subject: [PATCH] working on 3idx mo ints --- src/mo_basis/mos_complex.irp.f | 152 ++--------- src/mo_two_e_ints/EZFIO.cfg | 12 +- src/mo_two_e_ints/df_mo_ints.irp.f | 256 +++++++++++++++--- src/mo_two_e_ints/mo_bi_integrals.irp.f | 4 + src/utils_complex/MolPyscfToQPkpts.py | 182 +++---------- src/utils_complex/NEED | 1 + .../create_ezfio_complex_3idx.py | 28 +- 7 files changed, 313 insertions(+), 322 deletions(-) diff --git a/src/mo_basis/mos_complex.irp.f b/src/mo_basis/mos_complex.irp.f index 2e2f0786..e8c543ec 100644 --- a/src/mo_basis/mos_complex.irp.f +++ b/src/mo_basis/mos_complex.irp.f @@ -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) ] diff --git a/src/mo_two_e_ints/EZFIO.cfg b/src/mo_two_e_ints/EZFIO.cfg index fc1ff2e1..8cb039ea 100644 --- a/src/mo_two_e_ints/EZFIO.cfg +++ b/src/mo_two_e_ints/EZFIO.cfg @@ -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 diff --git a/src/mo_two_e_ints/df_mo_ints.irp.f b/src/mo_two_e_ints/df_mo_ints.irp.f index 5d85056b..7a0ab82a 100644 --- a/src/mo_two_e_ints/df_mo_ints.irp.f +++ b/src/mo_two_e_ints/df_mo_ints.irp.f @@ -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) diff --git a/src/mo_two_e_ints/mo_bi_integrals.irp.f b/src/mo_two_e_ints/mo_bi_integrals.irp.f index bdaa86c9..08fb82ba 100644 --- a/src/mo_two_e_ints/mo_bi_integrals.irp.f +++ b/src/mo_two_e_ints/mo_bi_integrals.irp.f @@ -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 diff --git a/src/utils_complex/MolPyscfToQPkpts.py b/src/utils_complex/MolPyscfToQPkpts.py index 55e4800d..5ebaa995 100644 --- a/src/utils_complex/MolPyscfToQPkpts.py +++ b/src/utils_complex/MolPyscfToQPkpts.py @@ -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,19 +430,17 @@ 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 - #idx2_cd = idx2_tri((c,d)) - for b, kb in enumerate(kpts): - if b > d: break - a = kconserv[b,c,d] - if a > d: continue - #if idx2_tri((a,b)) > idx2_cd: continue - #if ((c==d) and (a>b)): continue - ka = kpts[a] - with open(outfilename,'a') as outfile: + 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 a > d: continue + #if idx2_tri((a,b)) > idx2_cd: continue + #if ((c==d) and (a>b)): continue + ka = kpts[a] 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,21 +475,20 @@ 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 - #idx2_cd = idx2_tri((c,d)) - for b, kb in enumerate(kpts): - if b > d: break - a = kconserv[b,c,d] - if a > d: continue - #if idx2_tri((a,b)) > idx2_cd: continue - #if ((c==d) and (a>b)): continue - ka = kpts[a] + 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 a > d: continue + #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_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) diff --git a/src/utils_complex/NEED b/src/utils_complex/NEED index 173c6966..7b1c3363 100644 --- a/src/utils_complex/NEED +++ b/src/utils_complex/NEED @@ -1,2 +1,3 @@ ao_two_e_ints ao_one_e_ints +mo_two_e_ints diff --git a/src/utils_complex/create_ezfio_complex_3idx.py b/src/utils_complex/create_ezfio_complex_3idx.py index 31107376..9780b6e4 100755 --- a/src/utils_complex/create_ezfio_complex_3idx.py +++ b/src/utils_complex/create_ezfio_complex_3idx.py @@ -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,13 +136,24 @@ 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') -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') +# 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)) + 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