From a65eeab44e805bba1e07866b6ff93fd421eba789 Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Tue, 8 Mar 2022 17:48:25 -0600 Subject: [PATCH 01/28] starting 3-idx ints in qmcpack cholesky format --- src/ao_two_e_ints/EZFIO.cfg | 18 +++ src/ao_two_e_ints/cd_ao_ints.irp.f | 239 +++++++++++++++++++++++++++++ src/nuclei/EZFIO.cfg | 29 ++++ src/nuclei/kconserv_cplx.irp.f | 91 ++++++++++- 4 files changed, 375 insertions(+), 2 deletions(-) create mode 100644 src/ao_two_e_ints/cd_ao_ints.irp.f diff --git a/src/ao_two_e_ints/EZFIO.cfg b/src/ao_two_e_ints/EZFIO.cfg index 5b50f718..241d8a04 100644 --- a/src/ao_two_e_ints/EZFIO.cfg +++ b/src/ao_two_e_ints/EZFIO.cfg @@ -35,3 +35,21 @@ doc: Real part of the df integrals over AOs size: (2,ao_basis.ao_num_per_kpt,ao_basis.ao_num_per_kpt,ao_two_e_ints.df_num,nuclei.kpt_pair_num) interface: ezfio +[chol_num] +type: integer +doc: number of cholesky vecs for each kpt +size: (nuclei.unique_kpt_num) +interface: ezfio + +[chol_num_max] +type: integer +doc: max number of cholesky vecs +default: =maxval(ao_two_e_ints.chol_num) +interface: ezfio + +[chol_ao_integrals_complex] +type: double precision +doc: Cholesky decomposed integrals over AOs +size: (2,ao_basis.ao_num_per_kpt,ao_basis.ao_num_per_kpt,ao_two_e_ints.chol_num_max,nuclei.kpt_num,nuclei.unique_kpt_num) +interface: ezfio + diff --git a/src/ao_two_e_ints/cd_ao_ints.irp.f b/src/ao_two_e_ints/cd_ao_ints.irp.f new file mode 100644 index 00000000..7638aba9 --- /dev/null +++ b/src/ao_two_e_ints/cd_ao_ints.irp.f @@ -0,0 +1,239 @@ +BEGIN_PROVIDER [complex*16, chol_ao_integrals_complex, (ao_num_per_kpt,ao_num_per_kpt,chol_num_max,kpt_num,chol_unique_kpt_num)] + implicit none + BEGIN_DOC + ! CD AO integrals + ! first two dims are AOs x AOs + ! 3rd dim is chol_vec (pad with zeros to max size to avoid dealing with ragged array) + ! 4th dim is over all kpts + ! last dim is over "unique" kpts (one for each pair of additive inverses modulo G) + END_DOC + integer :: i,j,k,l + + if (read_chol_ao_integrals) then + call ezfio_get_ao_two_e_ints_chol_ao_integrals_complex(chol_ao_integrals_complex) + print *, 'CD AO integrals read from disk' + else + print*,'CD AO integrals must be provided',irp_here + stop -1 + endif + + if (write_chol_ao_integrals) then + call ezfio_set_ao_two_e_ints_chol_ao_integrals_complex(chol_ao_integrals_complex) + print *, 'CD AO integrals written to disk' + endif + +END_PROVIDER + + +subroutine ao_map_fill_from_chol + use map_module + implicit none + BEGIN_DOC + ! TODO: below is copy/paste of DF code as placeholder; modify for CD + ! fill ao bielec integral map using 3-index cd 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_ao,j_ao,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 :: ao_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 + + ao_num_kpt_2 = ao_num_per_kpt * ao_num_per_kpt + + size_buffer = min(ao_num_per_kpt*ao_num_per_kpt*ao_num_per_kpt,16000000) + print*, 'Providing the ao_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(ao_num_per_kpt,ao_num_per_kpt,chol_num_max)) + + 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_ao=1,ao_num_per_kpt + do j_ao=1,ao_num_per_kpt + do i_df=1,df_num + ints_jl(i_ao,j_ao,i_df) = dconjg(df_ao_integrals_complex(j_ao,i_ao,i_df,kjkl2)) + enddo + enddo + enddo + else + ints_jl = df_ao_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_ao, j_ao, 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, ao_num_per_kpt, ao_num_kpt_2, & + !$OMP chol_num_max, chol_num, chol_unique_kpt_num, chol_kpt_map, & + !$OMP kl,kj,kjkl2,ints_jl, & + !$OMP kconserv, df_ao_integrals_complex, ao_integrals_threshold, ao_integrals_map, ao_integrals_map_2) + + allocate( & + ints_ik(ao_num_per_kpt,ao_num_per_kpt,df_num), & + ints_ikjl(ao_num_per_kpt,ao_num_per_kpt,ao_num_per_kpt,ao_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_ao=1,ao_num_per_kpt + do j_ao=1,ao_num_per_kpt + do i_df=1,df_num + ints_ik(i_ao,j_ao,i_df) = dconjg(df_ao_integrals_complex(j_ao,i_ao,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_ao_integrals_complex(:,:,:,kikk2) + endif + + call zgemm('N','T', ao_num_kpt_2, ao_num_kpt_2, df_num, & + (1.d0,0.d0), ints_ik, ao_num_kpt_2, & + ints_jl, ao_num_kpt_2, & + (0.d0,0.d0), ints_ikjl, ao_num_kpt_2) + + n_integrals_1=0 + n_integrals_2=0 + do il=1,ao_num_per_kpt + l=il+(kl-1)*ao_num_per_kpt + do ij=1,ao_num_per_kpt + j=ij+(kj-1)*ao_num_per_kpt + if (j>l) exit + call idx2_tri_int(j,l,jl2) + do ik=1,ao_num_per_kpt + k=ik+(kk-1)*ao_num_per_kpt + if (k>l) exit + do ii=1,ao_num_per_kpt + i=ii+(ki-1)*ao_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) < ao_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 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 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 insert_into_ao_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1) + endif + if (n_integrals_2 > 0) then + 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(ao_integrals_map),'+',map_mb(ao_integrals_map_2),'MB' + endif + + enddo !kl + deallocate( ints_jl ) + + call map_sort(ao_integrals_map) + call map_unique(ao_integrals_map) + call map_sort(ao_integrals_map_2) + call map_unique(ao_integrals_map_2) + !call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_complex_1',ao_integrals_map) + !call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_complex_2',ao_integrals_map_2) + !call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read') + + call wall_time(wall_2) + call cpu_time(cpu_2) + + integer*8 :: get_ao_map_size, ao_map_size + ao_map_size = get_ao_map_size() + + print*,'AO integrals provided:' + print*,' Size of AO map ', map_mb(ao_integrals_map),'+',map_mb(ao_integrals_map_2),'MB' + print*,' Number of AO integrals: ', ao_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 ao_map_fill_from_df + diff --git a/src/nuclei/EZFIO.cfg b/src/nuclei/EZFIO.cfg index b4599b72..a9aae231 100644 --- a/src/nuclei/EZFIO.cfg +++ b/src/nuclei/EZFIO.cfg @@ -60,3 +60,32 @@ type: integer doc: array containing information about k-point symmetry size: (nuclei.kpt_num,nuclei.kpt_num,nuclei.kpt_num) interface: ezfio + +[kpt_pair_map] +type: integer +doc: mapping from pairs of kpts to total per electron +size: (nuclei.kpt_num,nuclei.kpt_num) +interface: ezfio + +[kpt_inv] +type: integer +doc: additive inverse for each kpt +size: (nuclei.kpt_num) +interface: ezfio + +[kpt_sparse_map] +type: integer +doc: mapping from kpt idx to unique idx, negative for conj. transp. +size: (nuclei.kpt_num) +interface: ezfio + +[unique_kpt_num] +type: integer +doc: number of pairs of kpts that are additive inverses (mod G) +interface: ezfio, provider + +[io_kpt_symnm] +doc: Read/Write kpt_symm arrays from/to disk [ Write | Read | None ] +type: Disk_access +interface: ezfio,provider,ocaml +default: None diff --git a/src/nuclei/kconserv_cplx.irp.f b/src/nuclei/kconserv_cplx.irp.f index 616ba779..a3da99bc 100644 --- a/src/nuclei/kconserv_cplx.irp.f +++ b/src/nuclei/kconserv_cplx.irp.f @@ -21,8 +21,9 @@ BEGIN_PROVIDER [integer, kconserv, (kpt_num,kpt_num,kpt_num)] call ezfio_get_nuclei_kconserv(kconserv) print *, 'kconserv read from disk' else - print*,'kconserv must be provided' - stop -1 + call set_kconserv(kconserv) + !print*,'kconserv must be provided' + !stop -1 endif if (write_kconserv) then call ezfio_set_nuclei_kconserv(kconserv) @@ -30,6 +31,76 @@ BEGIN_PROVIDER [integer, kconserv, (kpt_num,kpt_num,kpt_num)] endif END_PROVIDER +BEGIN_PROVIDER [integer, kpt_pair_map, (kpt_num,kpt_num)] + implicit none + BEGIN_DOC + ! Information about k-point symmetry + ! + ! for k-points I,K: kpt_pair_map(I,K) = \alpha + ! where Q_{\alpha} = k_I - k_K + ! + END_DOC + + if (read_kpt_symm) then + call ezfio_get_nuclei_kpt_pair_map(kpt_pair_map) + print *, 'kpt_pair_map read from disk' + else + print*,'kpt_pair_map must be provided' + stop -1 + endif + if (write_kpt_symm) then + call ezfio_set_nuclei_kpt_pair_map(kpt_pair_map) + print *, 'kpt_pair_map written to disk' + endif +END_PROVIDER + +BEGIN_PROVIDER [integer, kpt_inv, (kpt_num)] + implicit none + BEGIN_DOC + ! Information about k-point symmetry + ! + ! for k-point I: kpt_inv(I) = K + ! where k_I + k_K = 0 (mod G) + ! + END_DOC + + if (read_kpt_symm) then + call ezfio_get_nuclei_kpt_inv(kpt_inv) + print *, 'kpt_inv read from disk' + else + print*,'kpt_inv must be provided' + stop -1 + endif + if (write_kpt_symm) then + call ezfio_set_nuclei_kpt_inv(kpt_inv) + print *, 'kpt_inv written to disk' + endif +END_PROVIDER + +BEGIN_PROVIDER [integer, kpt_sparse_map, (kpt_num)] + implicit none + BEGIN_DOC + ! Information about k-point symmetry + ! + ! for k-point I: if kpt_sparse_map(I) = j + ! if j>0: data for k_I is stored at index j in chol_ints + ! if j<0: data for k_I is conj. transp. of data at index j in chol_{ao,mo}_integrals_complex + ! + END_DOC + + if (read_kpt_symm) then + call ezfio_get_nuclei_kpt_sparse_map(kpt_sparse_map) + print *, 'kpt_sparse_map read from disk' + else + print*,'kpt_sparse_map must be provided' + stop -1 + endif + if (write_kpt_symm) then + call ezfio_set_nuclei_kpt_sparse_map(kpt_sparse_map) + print *, 'kpt_sparse_map written to disk' + endif +END_PROVIDER + subroutine double_allowed_kpts(kh1,kh2,kp1,kp2,is_allowed) implicit none integer, intent(in) :: kh1,kh2,kp1,kp2 @@ -38,3 +109,19 @@ subroutine double_allowed_kpts(kh1,kh2,kp1,kp2,is_allowed) is_allowed = (kconserv(kh1,kh2,kp1) == kp2) end subroutine +subroutine set_kconserv(kcon) + implicit none + integer, intent(out) :: kcon(kpt_num,kpt_num,kpt_num) + integer :: i,j,k,qij + + do i=1,kpt_num + do k=1,kpt_num + ! Q = k_I - k_K + qik = kpt_pair_map(i,k) + do j=1,kpt_num + ! k_L = k_J - (-(k_I - k_K)) + kcon(i,j,k) = kpt_pair_map(j,kpt_inv(qik)) + enddo + enddo + enddo +end subroutine From 3fa67de8e62fb51f939f1eba44ac6d3da7b1043f Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Wed, 9 Mar 2022 13:35:15 -0600 Subject: [PATCH 02/28] renaming to match qmcpack --- src/nuclei/EZFIO.cfg | 4 ++-- src/nuclei/kconserv_cplx.irp.f | 44 +++++++++++++++++++++------------- 2 files changed, 29 insertions(+), 19 deletions(-) diff --git a/src/nuclei/EZFIO.cfg b/src/nuclei/EZFIO.cfg index a9aae231..bf5c562c 100644 --- a/src/nuclei/EZFIO.cfg +++ b/src/nuclei/EZFIO.cfg @@ -61,13 +61,13 @@ doc: array containing information about k-point symmetry size: (nuclei.kpt_num,nuclei.kpt_num,nuclei.kpt_num) interface: ezfio -[kpt_pair_map] +[qktok2] type: integer doc: mapping from pairs of kpts to total per electron size: (nuclei.kpt_num,nuclei.kpt_num) interface: ezfio -[kpt_inv] +[minusk] type: integer doc: additive inverse for each kpt size: (nuclei.kpt_num) diff --git a/src/nuclei/kconserv_cplx.irp.f b/src/nuclei/kconserv_cplx.irp.f index a3da99bc..e0958bec 100644 --- a/src/nuclei/kconserv_cplx.irp.f +++ b/src/nuclei/kconserv_cplx.irp.f @@ -31,49 +31,49 @@ BEGIN_PROVIDER [integer, kconserv, (kpt_num,kpt_num,kpt_num)] endif END_PROVIDER -BEGIN_PROVIDER [integer, kpt_pair_map, (kpt_num,kpt_num)] +BEGIN_PROVIDER [integer, qktok2, (kpt_num,kpt_num)] implicit none BEGIN_DOC ! Information about k-point symmetry ! - ! for k-points I,K: kpt_pair_map(I,K) = \alpha + ! for k-points I,K: qktok2(K,I) = \alpha ! where Q_{\alpha} = k_I - k_K ! END_DOC if (read_kpt_symm) then - call ezfio_get_nuclei_kpt_pair_map(kpt_pair_map) - print *, 'kpt_pair_map read from disk' + call ezfio_get_nuclei_qktok2(qktok2) + print *, 'qktok2 read from disk' else - print*,'kpt_pair_map must be provided' + print*,'qktok2 must be provided' stop -1 endif if (write_kpt_symm) then - call ezfio_set_nuclei_kpt_pair_map(kpt_pair_map) - print *, 'kpt_pair_map written to disk' + call ezfio_set_nuclei_qktok2(qktok2) + print *, 'qktok2 written to disk' endif END_PROVIDER -BEGIN_PROVIDER [integer, kpt_inv, (kpt_num)] +BEGIN_PROVIDER [integer, minusk, (kpt_num)] implicit none BEGIN_DOC ! Information about k-point symmetry ! - ! for k-point I: kpt_inv(I) = K + ! for k-point I: minusk(I) = K ! where k_I + k_K = 0 (mod G) ! END_DOC if (read_kpt_symm) then - call ezfio_get_nuclei_kpt_inv(kpt_inv) - print *, 'kpt_inv read from disk' + call ezfio_get_nuclei_minusk(minusk) + print *, 'minusk read from disk' else - print*,'kpt_inv must be provided' + print*,'minusk must be provided' stop -1 endif if (write_kpt_symm) then - call ezfio_set_nuclei_kpt_inv(kpt_inv) - print *, 'kpt_inv written to disk' + call ezfio_set_nuclei_minusk(minusk) + print *, 'minusk written to disk' endif END_PROVIDER @@ -85,6 +85,16 @@ BEGIN_PROVIDER [integer, kpt_sparse_map, (kpt_num)] ! for k-point I: if kpt_sparse_map(I) = j ! if j>0: data for k_I is stored at index j in chol_ints ! if j<0: data for k_I is conj. transp. of data at index j in chol_{ao,mo}_integrals_complex + ! + ! if we have h5 data stored under L[i]: + ! count=1 + ! do i=1,N_L + ! kpt_sparse_map(i)=count + ! if (minusk(i) != i) then + ! kpt_sparse_map(minusk(i)) = -count + ! endif + ! count += 1 + ! enddo ! END_DOC @@ -112,15 +122,15 @@ end subroutine subroutine set_kconserv(kcon) implicit none integer, intent(out) :: kcon(kpt_num,kpt_num,kpt_num) - integer :: i,j,k,qij + integer :: i,j,k,qik do i=1,kpt_num do k=1,kpt_num ! Q = k_I - k_K - qik = kpt_pair_map(i,k) + qik = qktok2(k,i) do j=1,kpt_num ! k_L = k_J - (-(k_I - k_K)) - kcon(i,j,k) = kpt_pair_map(j,kpt_inv(qik)) + kcon(i,j,k) = qktok2(minusk(j),qik) enddo enddo enddo From 9d4d03cdc68e5cc4e21b5e203a2115383fc21137 Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Wed, 9 Mar 2022 13:35:41 -0600 Subject: [PATCH 03/28] started ao map fill from cholesky need to check conj/transp correctness loops should be restructured ao map fill shouldn't actually be used, but might be good for testing, and mos will work the same way --- src/ao_two_e_ints/cd_ao_ints.irp.f | 62 +++++++++++++++++------------- 1 file changed, 36 insertions(+), 26 deletions(-) diff --git a/src/ao_two_e_ints/cd_ao_ints.irp.f b/src/ao_two_e_ints/cd_ao_ints.irp.f index 7638aba9..b9b602a4 100644 --- a/src/ao_two_e_ints/cd_ao_ints.irp.f +++ b/src/ao_two_e_ints/cd_ao_ints.irp.f @@ -37,7 +37,7 @@ subroutine ao_map_fill_from_chol integer :: ki,kk,kj,kl integer :: ii,ik,ij,il integer :: kikk2,kjkl2,jl2,ik2 - integer :: i_ao,j_ao,i_df + integer :: i_ao,j_ao,i_cd complex*16,allocatable :: ints_ik(:,:,:), ints_jl(:,:,:), ints_ikjl(:,:,:,:) @@ -59,7 +59,7 @@ subroutine ao_map_fill_from_chol ao_num_kpt_2 = ao_num_per_kpt * ao_num_per_kpt size_buffer = min(ao_num_per_kpt*ao_num_per_kpt*ao_num_per_kpt,16000000) - print*, 'Providing the ao_bielec integrals from 3-index df integrals' + print*, 'Providing the ao_bielec integrals from 3-index cholesky 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 @@ -70,34 +70,43 @@ subroutine ao_map_fill_from_chol allocate( ints_jl(ao_num_per_kpt,ao_num_per_kpt,chol_num_max)) wall_0 = wall_1 - do kl=1, kpt_num - do kj=1, kl + !TODO: change loops so that we only iterate over "correct" slices (i.e. ik block is stored directly, not as conj. transp.) + ! possible cases for (ik,jl) are (+,+), (+,-), (-,+), (-,-) + ! where + is the slice used as stored, and - is the conj. transp. of the stored data + ! (+,+) and (-,-) give the same information; we should always use (+,+) + ! (+,-) and (-,+) give the same information; we should always use (+,-) + do kQ = 1, kpt_num + do kl = 1, kpt_num + kj = qktok2(kQ,kl) + assert(kQ == qktok2(kj,kl)) + if (kj>kl) cycle call idx2_tri_int(kj,kl,kjkl2) - if (kj < kl) then + !TODO: verify the kj, kl as 4th index in expressions below + if (kpt_sparse_map(kQ) > 0) then + ints_jl = chol_ao_integrals_complex(:,:,:,kl,kpt_sparse_map(kQ)) + else do i_ao=1,ao_num_per_kpt do j_ao=1,ao_num_per_kpt - do i_df=1,df_num - ints_jl(i_ao,j_ao,i_df) = dconjg(df_ao_integrals_complex(j_ao,i_ao,i_df,kjkl2)) + do i_cd=1,chol_num_max + ints_jl(i_ao,j_ao,i_cd) = dconjg(chol_ao_integrals_complex(j_ao,i_ao,i_cd,kj,-kpt_sparse_map(kQ))) enddo enddo enddo - else - ints_jl = df_ao_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_ao, j_ao, i_df, & + !$OMP PARALLEL PRIVATE(i,k,j,l,ki,kk,ii,ik,ij,il,kikk2,jl2,ik2,kQ, & + !$OMP ints_ik, ints_ikjl, i_ao, j_ao, i_cd, & !$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, ao_num_per_kpt, ao_num_kpt_2, & - !$OMP chol_num_max, chol_num, chol_unique_kpt_num, chol_kpt_map, & + !$OMP chol_num_max, chol_num, chol_unique_kpt_num, kpt_sparse_map, qktok2, minusk, & !$OMP kl,kj,kjkl2,ints_jl, & - !$OMP kconserv, df_ao_integrals_complex, ao_integrals_threshold, ao_integrals_map, ao_integrals_map_2) + !$OMP kconserv, chol_ao_integrals_complex, ao_integrals_threshold, ao_integrals_map, ao_integrals_map_2) allocate( & - ints_ik(ao_num_per_kpt,ao_num_per_kpt,df_num), & + ints_ik(ao_num_per_kpt,ao_num_per_kpt,chol_num_max), & ints_ikjl(ao_num_per_kpt,ao_num_per_kpt,ao_num_per_kpt,ao_num_per_kpt), & buffer_i_1(size_buffer), & buffer_i_2(size_buffer), & @@ -107,25 +116,26 @@ subroutine ao_map_fill_from_chol !$OMP DO SCHEDULE(guided) do kk=1,kl - ki=kconserv(kl,kk,kj) + ki = qktok2(minusk(kk),kQ) + assert(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 + !TODO: check this! (ki, kk slice index and transpose/notranspose) + if (kpt_sparse_map(kQ) > 0) then + ints_ik = chol_ao_integrals_complex(:,:,:,ki,kpt_sparse_map(kQ)) + else do i_ao=1,ao_num_per_kpt do j_ao=1,ao_num_per_kpt - do i_df=1,df_num - ints_ik(i_ao,j_ao,i_df) = dconjg(df_ao_integrals_complex(j_ao,i_ao,i_df,kikk2)) + do i_cd=1,chol_num_max + ints_jl(i_ao,j_ao,i_cd) = dconjg(chol_ao_integrals_complex(j_ao,i_ao,i_cd,kk,-kpt_sparse_map(kQ))) 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_ao_integrals_complex(:,:,:,kikk2) endif - call zgemm('N','T', ao_num_kpt_2, ao_num_kpt_2, df_num, & + call zgemm('N','T', ao_num_kpt_2, ao_num_kpt_2, chol_num(kQ), & (1.d0,0.d0), ints_ik, ao_num_kpt_2, & ints_jl, ao_num_kpt_2, & (0.d0,0.d0), ints_ikjl, ao_num_kpt_2) @@ -204,15 +214,15 @@ subroutine ao_map_fill_from_chol buffer_values_2 & ) !$OMP END PARALLEL - enddo !kj + enddo !kl 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 ', & + print*, 100.*float(kQ)/float(kpt_num), '% in ', & wall_2-wall_1,'s',map_mb(ao_integrals_map),'+',map_mb(ao_integrals_map_2),'MB' endif - enddo !kl + enddo !kQ deallocate( ints_jl ) call map_sort(ao_integrals_map) @@ -235,5 +245,5 @@ subroutine ao_map_fill_from_chol 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 ao_map_fill_from_df +end subroutine ao_map_fill_from_chol From 1fb24382b570fa5b304c723ca8bc9e42ebbfe8e2 Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Wed, 9 Mar 2022 13:43:34 -0600 Subject: [PATCH 04/28] starting mo cholesky ints --- src/ao_two_e_ints/EZFIO.cfg | 6 ++++++ src/mo_two_e_ints/EZFIO.cfg | 12 ++++++++++++ src/mo_two_e_ints/mo_bi_integrals.irp.f | 6 ++++++ 3 files changed, 24 insertions(+) diff --git a/src/ao_two_e_ints/EZFIO.cfg b/src/ao_two_e_ints/EZFIO.cfg index 241d8a04..0dcb769e 100644 --- a/src/ao_two_e_ints/EZFIO.cfg +++ b/src/ao_two_e_ints/EZFIO.cfg @@ -47,6 +47,12 @@ doc: max number of cholesky vecs default: =maxval(ao_two_e_ints.chol_num) interface: ezfio +[io_chol_ao_integrals] +type: Disk_access +doc: Read/Write chol |AO| integrals from/to disk [ Write | Read | None ] +interface: ezfio,provider,ocaml +default: None + [chol_ao_integrals_complex] type: double precision doc: Cholesky decomposed integrals over AOs diff --git a/src/mo_two_e_ints/EZFIO.cfg b/src/mo_two_e_ints/EZFIO.cfg index c708792f..e8094c65 100644 --- a/src/mo_two_e_ints/EZFIO.cfg +++ b/src/mo_two_e_ints/EZFIO.cfg @@ -29,3 +29,15 @@ 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 +[io_chol_mo_integrals] +type: Disk_access +doc: Read/Write chol |MO| integrals from/to disk [ Write | Read | None ] +interface: ezfio,provider,ocaml +default: None + +[chol_mo_integrals_complex] +type: double precision +doc: Cholesky decomposed integrals over MOs +size: (2,mo_basis.mo_num_per_kpt,mo_basis.mo_num_per_kpt,ao_two_e_ints.chol_num_max,nuclei.kpt_num,nuclei.unique_kpt_num) +interface: ezfio + 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 e2447d6b..68af93bb 100644 --- a/src/mo_two_e_ints/mo_bi_integrals.irp.f +++ b/src/mo_two_e_ints/mo_bi_integrals.irp.f @@ -47,6 +47,12 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ] !call mo_map_fill_from_df_single call mo_map_fill_from_df_dot return + else if (read_chol_mo_integrals.or.read_chol_ao_integrals) then + PROVIDE chol_mo_integrals_complex + !call mo_map_fill_from_chol + !call mo_map_fill_from_chol_single + call mo_map_fill_from_chol_dot + return else PROVIDE ao_two_e_integrals_in_map endif From 30483cf0fa7081d1a0af7523a83f82129ad3bb5b Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Wed, 9 Mar 2022 14:30:36 -0600 Subject: [PATCH 05/28] cd ao mo transform --- src/mo_two_e_ints/cd_mo_ints.irp.f | 313 +++++++++++++++++++++++++++++ 1 file changed, 313 insertions(+) create mode 100644 src/mo_two_e_ints/cd_mo_ints.irp.f diff --git a/src/mo_two_e_ints/cd_mo_ints.irp.f b/src/mo_two_e_ints/cd_mo_ints.irp.f new file mode 100644 index 00000000..671edaba --- /dev/null +++ b/src/mo_two_e_ints/cd_mo_ints.irp.f @@ -0,0 +1,313 @@ +BEGIN_PROVIDER [complex*16, chol_mo_integrals_complex, (mo_num_per_kpt,mo_num_per_kpt,chol_num_max,kpt_num,chol_unique_kpt_num)] + implicit none + BEGIN_DOC + ! df MO integrals + END_DOC + integer :: i,j,k,l + + if (read_chol_mo_integrals) then + call ezfio_get_mo_two_e_ints_chol_mo_integrals_complex(chol_mo_integrals_complex) + print *, 'CD MO integrals read from disk' + else + call chol_mo_from_chol_ao(chol_mo_integrals_complex,chol_ao_integrals_complex,mo_num_per_kpt,ao_num_per_kpt, & + chol_num_max,kpt_num,chol_unique_kpt_num) + endif + + if (write_df_mo_integrals) then + 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_chol_dot + use map_module + implicit none + BEGIN_DOC + ! TODO: implement this (below is just copied from df as placeholder) + ! fill mo bielec integral map using 3-index cd integrals + END_DOC + + integer :: i,k,j,l,mu + 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(:,:,:) + + complex*16 :: integral,mjl,mik + 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 + !complex*16, external :: zdotc + complex*16, external :: zdotu + + 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(df_num,mo_num_per_kpt,mo_num_per_kpt)) + allocate( ints_ik(df_num,mo_num_per_kpt,mo_num_per_kpt)) + + 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_df,i_mo,j_mo) = dconjg(df_mo_integrals_complex(j_mo,i_mo,i_df,kjkl2)) + enddo + enddo + enddo + else + 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_df,i_mo,j_mo) = df_mo_integrals_complex(i_mo,j_mo,i_df,kjkl2) + enddo + enddo + enddo + endif + + do kk=1,kl + ki=kconserv(kl,kk,kj) + if (ki>kl) cycle + call idx2_tri_int(ki,kk,kikk2) + 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_df,i_mo,j_mo) = 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 + 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_df,i_mo,j_mo) = df_mo_integrals_complex(i_mo,j_mo,i_df,kikk2) + enddo + enddo + enddo + endif + + !$OMP PARALLEL PRIVATE(i,k,j,l,ii,ik,ij,il,jl2,ik2, & + !$OMP mu, mik, mjl, & + !$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 ki,kk,kikk2,ints_ik, & + !$OMP kconserv, df_mo_integrals_complex, mo_integrals_threshold, & + !$OMP mo_integrals_map, mo_integrals_map_2) + + allocate( & + buffer_i_1(size_buffer), & + buffer_i_2(size_buffer), & + buffer_values_1(size_buffer), & + buffer_values_2(size_buffer) & + ) + + n_integrals_1=0 + n_integrals_2=0 + !$OMP DO SCHEDULE(guided) + 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 = zdotc(df_num,ints_jl(1,ij,il),1,ints_ik(1,ii,ik),1) + integral = zdotu(df_num,ints_jl(1,ij,il),1,ints_ik(1,ii,ik),1) +! 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_mo_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1,mo_integrals_threshold) + 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_mo_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_2,mo_integrals_threshold) + n_integrals_2 = 0 + endif + endif + + enddo !ii + enddo !ik + enddo !ij + enddo !il + !$OMP END DO NOWAIT + + if (n_integrals_1 > 0) then + call map_append(mo_integrals_map, buffer_i_1, buffer_values_1, n_integrals_1) + !call insert_into_mo_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1,mo_integrals_threshold) + 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_mo_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_2,mo_integrals_threshold) + endif + deallocate( & + buffer_i_1, & + buffer_i_2, & + buffer_values_1, & + buffer_values_2 & + ) + !$OMP END PARALLEL + enddo !kk + 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,ints_ik ) + + 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_merge(mo_integrals_map) + !call map_merge(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_chol_dot + + +subroutine chol_mo_from_chol_ao(cd_mo,cd_ao,n_mo,n_ao,n_cd,n_k,n_unique_k) + use map_module + implicit none + BEGIN_DOC + ! create 3-idx mo ints from 3-idx ao ints + END_DOC + integer,intent(in) :: n_mo,n_ao,n_cd,n_k,n_unique_k + complex*16,intent(out) :: cd_mo(n_mo,n_mo,n_cd,n_k,n_unique_k) + complex*16,intent(in) :: cd_ao(n_ao,n_ao,n_cd,n_k,n_unique_k) + integer :: ki,kk,mu,kQ,Q_idx + complex*16,allocatable :: coef_i(:,:), coef_k(:,:), ints_ik(:,:), ints_tmp(:,:) + double precision :: wall_1,wall_2,cpu_1,cpu_2 + + print*,'providing 3-index MO integrals from 3-index AO integrals' + + cd_mo = 0.d0 + + call wall_time(wall_1) + call cpu_time(cpu_1) + allocate( & + coef_i(n_ao,n_mo),& + coef_k(n_ao,n_mo),& + ints_ik(n_ao,n_ao),& + ints_tmp(n_mo,n_ao)& + ) + + do ki=1, kpt_num + coef_i = mo_coef_complex_kpts(:,:,ki) + do kk=1, kpt_num + coef_k = mo_coef_complex_kpts(:,:,kk) + kQ = qktok2(kk,ki) + Q_idx = kpt_sparse_map(kQ) + if (Q_idx < 0) cycle + + do mu=1, chol_num(kQ) + ints_ik = cd_ao(:,:,mu,ki,Q_idx) + call zgemm('C','N',n_mo,n_ao,n_ao, & + (1.d0,0.d0), coef_i, n_ao, & + ints_ik, n_ao, & + (0.d0,0.d0), ints_tmp, n_mo) + + call zgemm('N','N',n_mo,n_mo,n_ao, & + (1.d0,0.d0), ints_tmp, n_mo, & + coef_k, n_ao, & + (0.d0,0.d0), cd_mo(:,:,mu,ki,Q_idx), n_mo) + enddo + enddo + call wall_time(wall_2) + print*,100.*float(ki)/kpt_num, '% in ', & + wall_2-wall_1, 's' + enddo + + deallocate( & + coef_i, & + coef_k, & + ints_ik, & + ints_tmp & + ) + call wall_time(wall_2) + call cpu_time(cpu_2) + print*,' 3-idx MO provided' + 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 chol_mo_from_chol_ao From 11b3dee6d075d654c51e33b17c1096cbe3afd235 Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Wed, 9 Mar 2022 14:48:09 -0600 Subject: [PATCH 06/28] initial implementation of MO map fill from cd needs testing (verify indexing/conj.transp.) --- src/mo_two_e_ints/cd_mo_ints.irp.f | 65 +++++++++++++++++------------- 1 file changed, 37 insertions(+), 28 deletions(-) diff --git a/src/mo_two_e_ints/cd_mo_ints.irp.f b/src/mo_two_e_ints/cd_mo_ints.irp.f index 671edaba..f270397d 100644 --- a/src/mo_two_e_ints/cd_mo_ints.irp.f +++ b/src/mo_two_e_ints/cd_mo_ints.irp.f @@ -1,7 +1,7 @@ BEGIN_PROVIDER [complex*16, chol_mo_integrals_complex, (mo_num_per_kpt,mo_num_per_kpt,chol_num_max,kpt_num,chol_unique_kpt_num)] implicit none BEGIN_DOC - ! df MO integrals + ! CD MO integrals END_DOC integer :: i,j,k,l @@ -13,9 +13,9 @@ BEGIN_PROVIDER [complex*16, chol_mo_integrals_complex, (mo_num_per_kpt,mo_num_pe chol_num_max,kpt_num,chol_unique_kpt_num) endif - if (write_df_mo_integrals) then - call ezfio_set_mo_two_e_ints_df_mo_integrals_complex(df_mo_integrals_complex) - print *, 'df MO integrals written to disk' + if (write_chol_mo_integrals) then + call ezfio_set_mo_two_e_ints_chol_mo_integrals_complex(chol_mo_integrals_complex) + print *, 'CD MO integrals written to disk' endif END_PROVIDER @@ -24,7 +24,7 @@ subroutine mo_map_fill_from_chol_dot use map_module implicit none BEGIN_DOC - ! TODO: implement this (below is just copied from df as placeholder) + ! TODO: verify correct indexing and conj.transp. ! fill mo bielec integral map using 3-index cd integrals END_DOC @@ -32,7 +32,8 @@ subroutine mo_map_fill_from_chol_dot integer :: ki,kk,kj,kl integer :: ii,ik,ij,il integer :: kikk2,kjkl2,jl2,ik2 - integer :: i_mo,j_mo,i_df + integer :: i_mo,j_mo,i_cd + integer :: kQ, Q_idx complex*16,allocatable :: ints_ik(:,:,:), ints_jl(:,:,:) @@ -56,7 +57,7 @@ subroutine mo_map_fill_from_chol_dot 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' + print*, 'Providing the mo_bielec integrals from 3-index CD 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 @@ -64,40 +65,47 @@ subroutine mo_map_fill_from_chol_dot call wall_time(wall_1) call cpu_time(cpu_1) - allocate( ints_jl(df_num,mo_num_per_kpt,mo_num_per_kpt)) - allocate( ints_ik(df_num,mo_num_per_kpt,mo_num_per_kpt)) + allocate( ints_jl(chol_num_max,mo_num_per_kpt,mo_num_per_kpt)) + allocate( ints_ik(chol_num_max,mo_num_per_kpt,mo_num_per_kpt)) wall_0 = wall_1 - do kl=1, kpt_num - do kj=1, kl + do kQ = 1, kpt_num + Q_idx = kpt_sparse_map(kQ) + do kl = 1, kpt_num + kj = qktok2(kQ,kl) + assert(kQ == qktok2(kj,kl)) + if (kj>kl) cycle call idx2_tri_int(kj,kl,kjkl2) - if (kj < kl) then + ints_jl = 0.d0 + if (Q_idx > 0) 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_df,i_mo,j_mo) = dconjg(df_mo_integrals_complex(j_mo,i_mo,i_df,kjkl2)) + do i_cd=1,chol_num(kQ) + ints_jl(i_cd,i_mo,j_mo) = chol_mo_integrals_complex(i_mo,j_mo,i_cd,kl,Q_idx) enddo enddo enddo else 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_df,i_mo,j_mo) = df_mo_integrals_complex(i_mo,j_mo,i_df,kjkl2) + do i_cd=1,chol_num(kQ) + ints_jl(i_cd,i_mo,j_mo) = dconjg(chol_mo_integrals_complex(j_mo,i_mo,i_cd,kj,-Q_idx)) enddo enddo enddo endif do kk=1,kl - ki=kconserv(kl,kk,kj) + ki = qktok2(minusk(kk),kQ) + assert(ki == kconserv(kl,kk,kj)) if (ki>kl) cycle call idx2_tri_int(ki,kk,kikk2) - if (ki < kk) then + ints_ik = 0.d0 + if (Q_idx > 0) 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_df,i_mo,j_mo) = dconjg(df_mo_integrals_complex(j_mo,i_mo,i_df,kikk2)) + do i_cd=1,chol_num(kQ) + ints_ik(i_cd,i_mo,j_mo) = chol_mo_integrals_complex(i_mo,j_mo,i_cd,ki,Q_idx) enddo enddo enddo @@ -105,8 +113,8 @@ subroutine mo_map_fill_from_chol_dot else 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_df,i_mo,j_mo) = df_mo_integrals_complex(i_mo,j_mo,i_df,kikk2) + do i_cd=1,chol_num(kQ) + ints_ik(i_cd,i_mo,j_mo) = dconjg(chol_mo_integrals_complex(j_mo,i_mo,i_cd,kk,-Q_idx)) enddo enddo enddo @@ -118,10 +126,11 @@ subroutine mo_map_fill_from_chol_dot !$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 SHARED(size_buffer, kpt_num, chol_num, mo_num_per_kpt, mo_num_kpt_2, & !$OMP kl,kj,kjkl2,ints_jl, & !$OMP ki,kk,kikk2,ints_ik, & - !$OMP kconserv, df_mo_integrals_complex, mo_integrals_threshold, & + !$OMP kQ, Q_idx, chol_num, & + !$OMP kconserv, chol_mo_integrals_complex, mo_integrals_threshold, & !$OMP mo_integrals_map, mo_integrals_map_2) allocate( & @@ -149,7 +158,7 @@ subroutine mo_map_fill_from_chol_dot call idx2_tri_int(i,k,ik2) if (ik2 > jl2) exit !integral = zdotc(df_num,ints_jl(1,ij,il),1,ints_ik(1,ii,ik),1) - integral = zdotu(df_num,ints_jl(1,ij,il),1,ints_ik(1,ii,ik),1) + integral = zdotu(chol_num(kQ),ints_jl(1,ij,il),1,ints_ik(1,ii,ik),1) ! print*,i,k,j,l,real(integral),imag(integral) if (cdabs(integral) < mo_integrals_threshold) then cycle @@ -209,15 +218,15 @@ subroutine mo_map_fill_from_chol_dot ) !$OMP END PARALLEL enddo !kk - enddo !kj + enddo !kl 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 ', & + print*, 100.*float(kQ)/float(kpt_num), '% in ', & wall_2-wall_1,'s',map_mb(mo_integrals_map),'+',map_mb(mo_integrals_map_2),'MB' endif - enddo !kl + enddo !kQ deallocate( ints_jl,ints_ik ) call map_sort(mo_integrals_map) From 682c5f54a3d615ad19791fbbd99d70d49dfd218b Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Wed, 9 Mar 2022 15:03:37 -0600 Subject: [PATCH 07/28] update comment --- src/ao_two_e_ints/cd_ao_ints.irp.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ao_two_e_ints/cd_ao_ints.irp.f b/src/ao_two_e_ints/cd_ao_ints.irp.f index b9b602a4..866117bf 100644 --- a/src/ao_two_e_ints/cd_ao_ints.irp.f +++ b/src/ao_two_e_ints/cd_ao_ints.irp.f @@ -29,7 +29,7 @@ subroutine ao_map_fill_from_chol use map_module implicit none BEGIN_DOC - ! TODO: below is copy/paste of DF code as placeholder; modify for CD + ! TODO: check indexing/conj.transp. of slices; restructure loops ! fill ao bielec integral map using 3-index cd integrals END_DOC From 7ce27f6ce83e1a0c55f5cfdf24651f450ee02deb Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Wed, 9 Mar 2022 16:31:13 -0600 Subject: [PATCH 08/28] fixed typo --- src/nuclei/EZFIO.cfg | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/nuclei/EZFIO.cfg b/src/nuclei/EZFIO.cfg index bf5c562c..79774ea5 100644 --- a/src/nuclei/EZFIO.cfg +++ b/src/nuclei/EZFIO.cfg @@ -84,7 +84,7 @@ type: integer doc: number of pairs of kpts that are additive inverses (mod G) interface: ezfio, provider -[io_kpt_symnm] +[io_kpt_symm] doc: Read/Write kpt_symm arrays from/to disk [ Write | Read | None ] type: Disk_access interface: ezfio,provider,ocaml From 618df9eed57dbc50bd90d27e33557667148e4a96 Mon Sep 17 00:00:00 2001 From: amandadumi Date: Tue, 3 May 2022 12:23:51 -0400 Subject: [PATCH 09/28] starting separate function with dummy atom --- bin/qp_convert_h5_to_ezfio | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/bin/qp_convert_h5_to_ezfio b/bin/qp_convert_h5_to_ezfio index 274c6131..fb09097b 100755 --- a/bin/qp_convert_h5_to_ezfio +++ b/bin/qp_convert_h5_to_ezfio @@ -9,6 +9,7 @@ Options: -o --output=EZFIO_DIR Produced directory by default is FILE.ezfio --noqmc don't include basis, cell, etc. for QMCPACK + -cd h5 contains cholesky decomposition informatin, these h5 result from RMG and the pyscf AFQMC converter of QMCPACK. """ from ezfio import ezfio @@ -138,6 +139,17 @@ def convert_mol(filename,qph5path): return +def convert_kpts_cd(filename,qph5path,qmcpack=True): + + ezfio.set_file(filename) + ezfio.set_nuclei_is_complex(True) + # Dummy atom since AFQMC h5 has no atom information + nucl_num = 1 + ezfio.set_nuclei_nucl_num(nucl_num) + ezfio.set_nuclei_nucl_charge([0.]*nucl_num) + ezfio.set_nuclei_nucl_coord( [ [0.], [0.], [0.] ]*nucl_num ) + ezfio.set_nuclei_nucl_label( ['He'] * nucl_num ) + def convert_kpts(filename,qph5path,qmcpack=True): ezfio.set_file(filename) ezfio.set_nuclei_is_complex(True) @@ -640,10 +652,16 @@ if __name__ == '__main__': if ARGUMENTS["--noqmc"]: qmcpack = False + if ARGUMENTS["--rmg"]: + rmg = True with h5py.File(FILE,'r') as qph5: do_kpts = ('kconserv' in qph5['nuclei'].keys()) if (do_kpts): print("converting HDF5 to EZFIO for periodic system") + if cd: + print("Using RMG and AFQMC h5") + convert_kpts_cd(EZFIO_FILE,FILE,qmcpack) + else: convert_kpts(EZFIO_FILE,FILE,qmcpack) else: print("converting HDF5 to EZFIO for molecular system") From 3a7e30af73b3dbcf62b46bb8f2cde586e6e2ef61 Mon Sep 17 00:00:00 2001 From: amandadumi Date: Wed, 4 May 2022 15:14:04 -0400 Subject: [PATCH 10/28] adding chol_num info and function framework --- .gitignore | 1 - bin/qp_convert_h5_to_ezfio | 130 ++++++++++++++++++++++++++++++++++++- 2 files changed, 129 insertions(+), 2 deletions(-) diff --git a/.gitignore b/.gitignore index 096e385b..7b248722 100644 --- a/.gitignore +++ b/.gitignore @@ -3,7 +3,6 @@ quantum_package_static.tar.gz build.ninja .ninja_log .ninja_deps -bin/ lib/ config/qp_create_ninja.pickle src/*/.gitignore diff --git a/bin/qp_convert_h5_to_ezfio b/bin/qp_convert_h5_to_ezfio index fb09097b..adfcbdf7 100755 --- a/bin/qp_convert_h5_to_ezfio +++ b/bin/qp_convert_h5_to_ezfio @@ -149,6 +149,134 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True): ezfio.set_nuclei_nucl_charge([0.]*nucl_num) ezfio.set_nuclei_nucl_coord( [ [0.], [0.], [0.] ]*nucl_num ) ezfio.set_nuclei_nucl_label( ['He'] * nucl_num ) + with h5py.File(qph5path,'r') as qph5: + kpt_num = qph5['Hamiltonian']['KPoints'].shape[0] + elec_alpha_num = qph5['Hamiltonian']['dims'][4] + elec_beta_num = qph5['Hamiltonian']['dims'][5] + orb_num = qph5['Hamiltonian']['dims'][3] + try: + is_ao = json.loads(qph5['metadata'][()].decode("utf-8").replace("'",'"'))['ortho_ao'] + if is_ao: + ao_num = orb_num + elif is_ao ==False: + mo_num = orb_num + else: + raise ValueError('Problem with ortho_ao key in metadata') + except: + raise UnicodeDecodeError('metadata not correctly parsed from HDF5 file') + + 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) + # don't multiply nuclei by kpt_num + # work in k-space, not in equivalent supercell + nucl_num_per_kpt = nucl_num + ezfio.set_nuclei_nucl_num(nucl_num_per_kpt) + # these are totals (kpt_num * num_per_kpt) + # need to change if we want to truncate orbital space within pyscf + ezfio.set_ao_basis_ao_num(ao_num) + ezfio.set_mo_basis_mo_num(mo_num) + ezfio.set_ao_basis_ao_num_per_kpt(ao_num//kpt_num) + ezfio.set_mo_basis_mo_num_per_kpt(mo_num//kpt_num) + ezfio.electrons_elec_alpha_num = elec_alpha_num + ezfio.electrons_elec_beta_num = elec_beta_num + ########################################## + # # + # Basis # + # # + ########################################## + #TODO + ########################################## + # # + # MOCoeff # + # # + ########################################## + #TODO + ########################################## + # # + # Integrals Mono # + # # + ########################################## + """ + with h5py.File(qph5path,'r') as qph5: + if is_ao: + kin_ao_reim= + ovlp_ao_reim= + ne_ao_reim= + + ezfio.set_ao_one_e_ints_ao_integrals_kinetic_kpts(kin_ao_reim) + ezfio.set_ao_one_e_ints_ao_integrals_overlap_kpts(ovlp_ao_reim) + ezfio.set_ao_one_e_ints_ao_integrals_n_e_kpts(ne_ao_reim) + + 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') + else: + kin_mo_reim= + ovlp_mo_reim= + ne_mo_reim= + + ezfio.set_mo_one_e_ints_mo_integrals_kinetic_kpts(kin_mo_reim) + ezfio.set_mo_one_e_ints_mo_integrals_overlap_kpts(ovlp_mo_reim) + ezfio.set_mo_one_e_ints_mo_integrals_n_e_kpts(ne_mo_reim) + + ezfio.set_mo_one_e_ints_io_mo_integrals_kinetic('Read') + ezfio.set_mo_one_e_ints_io_mo_integrals_overlap('Read') + ezfio.set_mo_one_e_ints_io_mo_integrals_n_e('Read') + """ + ########################################## + # # + # k-points # + # # + ########################################## + #TODO + with h5py.File(qph5path,'r') as qph5: + kconserv = qph5['nuclei/kconserv'][()].tolist() + + ezfio.set_nuclei_kconserv(kconserv) + ezfio.set_nuclei_io_kconserv('Read') + # qktok2 + # minuxk + # kpt_sparse_map + # unique_kpt_num + # io_kpt_symm + ########################################## + # # + # Integrals Bi # + # # + ########################################## + + # should this be in ao_basis? ao_two_e_ints? + with h5py.File(qph5path,'r') as qph5: + qph5[] + nchol_per_kpt = qph5['Hamiltonian']['NCholPerKP'][:] + nchol_per_kpt = nchol_per_kpt[nchol_per_kpt != 0] + ezfio.set_chol_num(nchol_per_kpt) + ezfio.set_chol_num_max(max(nchol_per_kpt)) + if is_ao: + ezfio.set_io_chol_ao_integrals('Read') + """ + df_num = qph5['ao_two_e_ints'].attrs['df_num'] + ezfio.set_ao_two_e_ints_df_num(df_num) + if 'df_ao_integrals' in qph5['ao_two_e_ints'].keys(): + dfao_reim=qph5['ao_two_e_ints/df_ao_integrals'][()].tolist() + ezfio.set_ao_two_e_ints_df_ao_integrals_complex(dfao_reim) + ezfio.set_ao_two_e_ints_io_df_ao_integrals('Read') + """ + else: + ezfio.set_io_chol_mo_integrals('Read') + """ + ezfio.set_io_chol_mo_integrals('Read') + df_num = qph5['ao_two_e_ints'].attrs['df_num'] + ezfio.set_ao_two_e_ints_df_num(df_num) + dfmo_reim=qph5['mo_two_e_ints/df_mo_integrals'][()].tolist() + ezfio.set_mo_two_e_ints_df_mo_integrals_complex(dfmo_reim) + ezfio.set_mo_two_e_ints_io_df_mo_integrals('Read') + """ + return + + + def convert_kpts(filename,qph5path,qmcpack=True): ezfio.set_file(filename) @@ -662,7 +790,7 @@ if __name__ == '__main__': print("Using RMG and AFQMC h5") convert_kpts_cd(EZFIO_FILE,FILE,qmcpack) else: - convert_kpts(EZFIO_FILE,FILE,qmcpack) + convert_kpts(EZFIO_FILE,FILE,qmcpack) else: print("converting HDF5 to EZFIO for molecular system") convert_mol(EZFIO_FILE,FILE) From 3f924903abf40b6b7cc06d7dddeaee373dae61f9 Mon Sep 17 00:00:00 2001 From: Shiv Upadhyay Date: Wed, 4 May 2022 16:03:11 -0400 Subject: [PATCH 11/28] Start ao chol ints --- bin/qp_convert_h5_to_ezfio | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/bin/qp_convert_h5_to_ezfio b/bin/qp_convert_h5_to_ezfio index adfcbdf7..0b0c80e8 100755 --- a/bin/qp_convert_h5_to_ezfio +++ b/bin/qp_convert_h5_to_ezfio @@ -251,10 +251,26 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True): qph5[] nchol_per_kpt = qph5['Hamiltonian']['NCholPerKP'][:] nchol_per_kpt = nchol_per_kpt[nchol_per_kpt != 0] + nchol_per_kpt_max = max(nchol_per_kpt) ezfio.set_chol_num(nchol_per_kpt) - ezfio.set_chol_num_max(max(nchol_per_kpt)) + ezfio.set_chol_num_max(nchol_per_kpt_max) if is_ao: + ao_num_per_kpt = ao_num//kpt_num ezfio.set_io_chol_ao_integrals('Read') + ao_chol_two_e_ints = np.zeros((2, ao_num_per_kpt, ao_num_per_kpt, nchol_per_kpt_max, kpt_num, len(nchol_per_kpt))) + + for i in len(nchol_per_kpt): + L = qph5['Hamiltonian']['KPFactorized'][f'L{i}'][:] + L.reshape(kpt_num, ao_num_per_kpt, ao_num_per_kpt, nchol_per_kpt[i], 2) + #(6, 5184, 2) + for cmplx in range(2): + for ao_idx_i in range(ao_num_per_kpt): + for ao_idx_j in range(ao_num_per_kpt): + for chol_idx in range(nchol_per_kpt[i]): + for kpt_idx in range(kpt_num): + ao_chol_two_e_ints[cmplx][ao_idx_i][ao_idx_j][chol_idx][kpt_idx][i] = L[kpt_idx][ao_idx_i][ao_idx_i][chol_idx][cmplx] + + #(2,ao_basis.ao_num_per_kpt,ao_basis.ao_num_per_kpt,ao_two_e_ints.chol_num_max,nuclei.kpt_num,nuclei.unique_kpt_num#) """ df_num = qph5['ao_two_e_ints'].attrs['df_num'] ezfio.set_ao_two_e_ints_df_num(df_num) From 4893703e9c4e17dad4cfb3bb5533b29f247e9546 Mon Sep 17 00:00:00 2001 From: Shiv Upadhyay Date: Wed, 4 May 2022 16:49:05 -0400 Subject: [PATCH 12/28] Reordering indicies AO chol integrals --- bin/qp_convert_h5_to_ezfio | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/bin/qp_convert_h5_to_ezfio b/bin/qp_convert_h5_to_ezfio index 0b0c80e8..ddda2acf 100755 --- a/bin/qp_convert_h5_to_ezfio +++ b/bin/qp_convert_h5_to_ezfio @@ -257,18 +257,26 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True): if is_ao: ao_num_per_kpt = ao_num//kpt_num ezfio.set_io_chol_ao_integrals('Read') - ao_chol_two_e_ints = np.zeros((2, ao_num_per_kpt, ao_num_per_kpt, nchol_per_kpt_max, kpt_num, len(nchol_per_kpt))) - + #ao_chol_two_e_ints = np.zeros((2, ao_num_per_kpt, ao_num_per_kpt, nchol_per_kpt_max, kpt_num, len(nchol_per_kpt))) + L_list = [] for i in len(nchol_per_kpt): L = qph5['Hamiltonian']['KPFactorized'][f'L{i}'][:] L.reshape(kpt_num, ao_num_per_kpt, ao_num_per_kpt, nchol_per_kpt[i], 2) + L = np.einsum("ijklm->ilkjm", A, B) + L_list.append(L) + #(6, 5184, 2) + """ for cmplx in range(2): for ao_idx_i in range(ao_num_per_kpt): for ao_idx_j in range(ao_num_per_kpt): for chol_idx in range(nchol_per_kpt[i]): for kpt_idx in range(kpt_num): - ao_chol_two_e_ints[cmplx][ao_idx_i][ao_idx_j][chol_idx][kpt_idx][i] = L[kpt_idx][ao_idx_i][ao_idx_i][chol_idx][cmplx] + ao_chol_two_e_ints[cmplx][ao_idx_i][ao_idx_j][chol_idx][kpt_idx][i] = L[kpt_idx][ao_idx_i][ao_idx_j][chol_idx][cmplx] + """ + ao_chol_two_e_ints = np.vstack(L_list) + ao_chol_two_e_ints = ao_chol_two_e_ints.transpose() + ezfio.set_chol_ao_integrals_complex(ao_chol_two_e_ints) #(2,ao_basis.ao_num_per_kpt,ao_basis.ao_num_per_kpt,ao_two_e_ints.chol_num_max,nuclei.kpt_num,nuclei.unique_kpt_num#) """ From e63c95c9edac5f59d93f9c3f00236873468422a5 Mon Sep 17 00:00:00 2001 From: amandadumi Date: Wed, 11 May 2022 15:35:44 -0400 Subject: [PATCH 13/28] adding mo 2e ints --- bin/qp_convert_h5_to_ezfio | 277 ++++++++++++++++++++----------------- 1 file changed, 148 insertions(+), 129 deletions(-) diff --git a/bin/qp_convert_h5_to_ezfio b/bin/qp_convert_h5_to_ezfio index ddda2acf..f1782867 100755 --- a/bin/qp_convert_h5_to_ezfio +++ b/bin/qp_convert_h5_to_ezfio @@ -30,23 +30,23 @@ def get_full_path(file_path): def convert_mol(filename,qph5path): ezfio.set_file(filename) ezfio.set_nuclei_is_complex(False) - + with h5py.File(qph5path,'r') as qph5: nucl_num = qph5['nuclei'].attrs['nucl_num'] ao_num = qph5['ao_basis'].attrs['ao_num'] mo_num = qph5['mo_basis'].attrs['mo_num'] elec_alpha_num = qph5['electrons'].attrs['elec_alpha_num'] elec_beta_num = qph5['electrons'].attrs['elec_beta_num'] - + ezfio.set_nuclei_nucl_num(nucl_num) - + ezfio.set_ao_basis_ao_num(ao_num) ezfio.set_mo_basis_mo_num(mo_num) ezfio.electrons_elec_alpha_num = elec_alpha_num ezfio.electrons_elec_beta_num = elec_beta_num - - - + + + ##ao_num = mo_num ##Important ! #import math @@ -56,42 +56,42 @@ def convert_mol(filename,qph5path): # #ezfio.electrons_elec_alpha_num = int(nelec_alpha_per_kpt * n_kpts) #ezfio.electrons_elec_beta_num = int(nelec_beta_per_kpt * n_kpts) - + #ezfio.electrons_elec_alpha_num = int(math.ceil(num_elec / 2.)) #ezfio.electrons_elec_beta_num = int(math.floor(num_elec / 2.)) - + #ezfio.set_utils_num_kpts(n_kpts) #ezfio.set_integrals_bielec_df_num(n_aux) - + #(old)Important #ezfio.set_nuclei_nucl_num(nucl_num) #ezfio.set_nuclei_nucl_charge([0.]*nucl_num) #ezfio.set_nuclei_nucl_coord( [ [0.], [0.], [0.] ]*nucl_num ) #ezfio.set_nuclei_nucl_label( ['He'] * nucl_num ) - - + + with h5py.File(qph5path,'r') as qph5: nucl_charge=qph5['nuclei/nucl_charge'][()].tolist() nucl_coord=qph5['nuclei/nucl_coord'][()].T.tolist() nucl_label=qph5['nuclei/nucl_label'][()].tolist() nuclear_repulsion = qph5['nuclei'].attrs['nuclear_repulsion'] - + ezfio.set_nuclei_nucl_charge(nucl_charge) ezfio.set_nuclei_nucl_coord(nucl_coord) if isinstance(nucl_label[0],bytes): nucl_label = list(map(lambda x:x.decode(),nucl_label)) ezfio.set_nuclei_nucl_label(nucl_label) - + ezfio.set_nuclei_io_nuclear_repulsion('Read') ezfio.set_nuclei_nuclear_repulsion(nuclear_repulsion) - - + + ########################################## # # # Basis # # # ########################################## - + with h5py.File(qph5path,'r') as qph5: do_pseudo = qph5['pseudo'].attrs['do_pseudo'] ezfio.set_pseudo_do_pseudo(do_pseudo) @@ -106,13 +106,13 @@ def convert_mol(filename,qph5path): ezfio.set_pseudo_pseudo_v_kl(qph5['pseudo/pseudo_v_kl'][()].tolist()) ezfio.set_pseudo_pseudo_dz_k(qph5['pseudo/pseudo_dz_k'][()].tolist()) ezfio.set_pseudo_pseudo_dz_kl(qph5['pseudo/pseudo_dz_kl'][()].tolist()) - + ########################################## # # # Basis # # # ########################################## - + with h5py.File(qph5path,'r') as qph5: #coeftmp = qph5['ao_basis/ao_coef'][()] #expotmp = qph5['ao_basis/ao_expo'][()] @@ -122,25 +122,25 @@ def convert_mol(filename,qph5path): ezfio.set_ao_basis_ao_power(qph5['ao_basis/ao_power'][()].tolist()) ezfio.set_ao_basis_ao_coef(qph5['ao_basis/ao_coef'][()].tolist()) ezfio.set_ao_basis_ao_expo(qph5['ao_basis/ao_expo'][()].tolist()) - - + + ########################################## # # # MO Coef # # # ########################################## - - + + with h5py.File(qph5path,'r') as qph5: mo_coef = qph5['mo_basis/mo_coef'][()].tolist() ezfio.set_mo_basis_mo_coef(mo_coef) #maybe fix qp so we don't need this? #ezfio.set_mo_basis_mo_coef([[i for i in range(mo_num)] * ao_num]) - + return def convert_kpts_cd(filename,qph5path,qmcpack=True): - + ezfio.set_file(filename) ezfio.set_nuclei_is_complex(True) # Dummy atom since AFQMC h5 has no atom information @@ -170,7 +170,7 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True): ezfio.set_nuclei_kpt_pair_num(kpt_pair_num) # don't multiply nuclei by kpt_num # work in k-space, not in equivalent supercell - nucl_num_per_kpt = nucl_num + nucl_num_per_kpt = nucl_num ezfio.set_nuclei_nucl_num(nucl_num_per_kpt) # these are totals (kpt_num * num_per_kpt) # need to change if we want to truncate orbital space within pyscf @@ -203,23 +203,23 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True): kin_ao_reim= ovlp_ao_reim= ne_ao_reim= - + ezfio.set_ao_one_e_ints_ao_integrals_kinetic_kpts(kin_ao_reim) ezfio.set_ao_one_e_ints_ao_integrals_overlap_kpts(ovlp_ao_reim) ezfio.set_ao_one_e_ints_ao_integrals_n_e_kpts(ne_ao_reim) - + 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') - else: + else: kin_mo_reim= ovlp_mo_reim= ne_mo_reim= - + ezfio.set_mo_one_e_ints_mo_integrals_kinetic_kpts(kin_mo_reim) ezfio.set_mo_one_e_ints_mo_integrals_overlap_kpts(ovlp_mo_reim) ezfio.set_mo_one_e_ints_mo_integrals_n_e_kpts(ne_mo_reim) - + ezfio.set_mo_one_e_ints_io_mo_integrals_kinetic('Read') ezfio.set_mo_one_e_ints_io_mo_integrals_overlap('Read') ezfio.set_mo_one_e_ints_io_mo_integrals_n_e('Read') @@ -229,26 +229,25 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True): # k-points # # # ########################################## - #TODO + #TODO with h5py.File(qph5path,'r') as qph5: kconserv = qph5['nuclei/kconserv'][()].tolist() - + ezfio.set_nuclei_kconserv(kconserv) ezfio.set_nuclei_io_kconserv('Read') # qktok2 # minuxk # kpt_sparse_map # unique_kpt_num - # io_kpt_symm + # io_kpt_symm ########################################## # # # Integrals Bi # # # ########################################## - + # should this be in ao_basis? ao_two_e_ints? with h5py.File(qph5path,'r') as qph5: - qph5[] nchol_per_kpt = qph5['Hamiltonian']['NCholPerKP'][:] nchol_per_kpt = nchol_per_kpt[nchol_per_kpt != 0] nchol_per_kpt_max = max(nchol_per_kpt) @@ -288,7 +287,6 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True): ezfio.set_ao_two_e_ints_io_df_ao_integrals('Read') """ else: - ezfio.set_io_chol_mo_integrals('Read') """ ezfio.set_io_chol_mo_integrals('Read') df_num = qph5['ao_two_e_ints'].attrs['df_num'] @@ -296,16 +294,38 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True): dfmo_reim=qph5['mo_two_e_ints/df_mo_integrals'][()].tolist() ezfio.set_mo_two_e_ints_df_mo_integrals_complex(dfmo_reim) ezfio.set_mo_two_e_ints_io_df_mo_integrals('Read') - """ - return + """ + mo_num_per_kpt = ao_num//kpt_num + ezfio.set_io_chol_mo_integrals('Read') + #ao_chol_two_e_ints = np.zeros((2, ao_num_per_kpt, ao_num_per_kpt, nchol_per_kpt_max, kpt_num, len(nchol_per_kpt))) + L_list = [] + for i in len(nchol_per_kpt): + L = qph5['Hamiltonian']['KPFactorized'][f'L{i}'][:] + L.reshape(kpt_num, mo_num_per_kpt, mo_num_per_kpt, nchol_per_kpt[i], 2) + L = np.einsum("ijklm->ilkjm", A, B) + L_list.append(L) + + #(6, 5184, 2) + """ + for cmplx in range(2): + for ao_idx_i in range(ao_num_per_kpt): + for ao_idx_j in range(ao_num_per_kpt): + for chol_idx in range(nchol_per_kpt[i]): + for kpt_idx in range(kpt_num): + ao_chol_two_e_ints[cmplx][ao_idx_i][ao_idx_j][chol_idx][kpt_idx][i] = L[kpt_idx][ao_idx_i][ao_idx_j][chol_idx][cmplx] + """ + mo_chol_two_e_ints = np.vstack(L_list) + mo_chol_two_e_ints = mo_chol_two_e_ints.transpose() + ezfio.set_chol_mo_integrals_complex(mo_chol_two_e_ints) + return + - def convert_kpts(filename,qph5path,qmcpack=True): ezfio.set_file(filename) ezfio.set_nuclei_is_complex(True) - + with h5py.File(qph5path,'r') as qph5: kpt_num = qph5['nuclei'].attrs['kpt_num'] nucl_num = qph5['nuclei'].attrs['nucl_num'] @@ -313,16 +333,16 @@ def convert_kpts(filename,qph5path,qmcpack=True): mo_num = qph5['mo_basis'].attrs['mo_num'] elec_alpha_num = qph5['electrons'].attrs['elec_alpha_num'] elec_beta_num = qph5['electrons'].attrs['elec_beta_num'] - + 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) - + # don't multiply nuclei by kpt_num # work in k-space, not in equivalent supercell - nucl_num_per_kpt = nucl_num + nucl_num_per_kpt = nucl_num ezfio.set_nuclei_nucl_num(nucl_num_per_kpt) - + # these are totals (kpt_num * num_per_kpt) # need to change if we want to truncate orbital space within pyscf ezfio.set_ao_basis_ao_num(ao_num) @@ -331,9 +351,9 @@ def convert_kpts(filename,qph5path,qmcpack=True): ezfio.set_mo_basis_mo_num_per_kpt(mo_num//kpt_num) ezfio.electrons_elec_alpha_num = elec_alpha_num ezfio.electrons_elec_beta_num = elec_beta_num - - - + + + ##ao_num = mo_num ##Important ! #import math @@ -343,57 +363,57 @@ def convert_kpts(filename,qph5path,qmcpack=True): # #ezfio.electrons_elec_alpha_num = int(nelec_alpha_per_kpt * n_kpts) #ezfio.electrons_elec_beta_num = int(nelec_beta_per_kpt * n_kpts) - + #ezfio.electrons_elec_alpha_num = int(math.ceil(num_elec / 2.)) #ezfio.electrons_elec_beta_num = int(math.floor(num_elec / 2.)) - + #ezfio.set_utils_num_kpts(n_kpts) #ezfio.set_integrals_bielec_df_num(n_aux) - + #(old)Important #ezfio.set_nuclei_nucl_num(nucl_num) #ezfio.set_nuclei_nucl_charge([0.]*nucl_num) #ezfio.set_nuclei_nucl_coord( [ [0.], [0.], [0.] ]*nucl_num ) #ezfio.set_nuclei_nucl_label( ['He'] * nucl_num ) - - + + with h5py.File(qph5path,'r') as qph5: nucl_charge=qph5['nuclei/nucl_charge'][()].tolist() nucl_coord=qph5['nuclei/nucl_coord'][()].T.tolist() nucl_label=qph5['nuclei/nucl_label'][()].tolist() nuclear_repulsion = qph5['nuclei'].attrs['nuclear_repulsion'] - + ezfio.set_nuclei_nucl_charge(nucl_charge) ezfio.set_nuclei_nucl_coord(nucl_coord) if isinstance(nucl_label[0],bytes): nucl_label = list(map(lambda x:x.decode(),nucl_label)) ezfio.set_nuclei_nucl_label(nucl_label) - + ezfio.set_nuclei_io_nuclear_repulsion('Read') ezfio.set_nuclei_nuclear_repulsion(nuclear_repulsion) - - + + ########################################## # # # Basis(Dummy) # # # ########################################## - + with h5py.File(qph5path,'r') as qph5: ezfio.set_ao_basis_ao_basis(qph5['ao_basis'].attrs['ao_basis']) ezfio.set_ao_basis_ao_nucl(qph5['ao_basis/ao_nucl'][()].tolist()) - - + + #Just need one (can clean this up later) ao_prim_num_max = 5 - + d = [ [0] *ao_prim_num_max]*ao_num ezfio.set_ao_basis_ao_prim_num([ao_prim_num_max]*ao_num) ezfio.set_ao_basis_ao_power(d) ezfio.set_ao_basis_ao_coef(d) ezfio.set_ao_basis_ao_expo(d) - - + + ########################################### ## # ## Pseudo # @@ -431,7 +451,7 @@ def convert_kpts(filename,qph5path,qmcpack=True): # ezfio.set_ao_basis_ao_coef(qph5['ao_basis/ao_coef'][()].tolist()) # ezfio.set_ao_basis_ao_expo(qph5['ao_basis/ao_expo'][()].tolist()) - + ########################################## # # # Basis(QMC) # @@ -448,7 +468,7 @@ def convert_kpts(filename,qph5path,qmcpack=True): ezfio.set_qmcpack_qmc_lbas(qph5['qmcpack/qmc_lbas'][()].tolist()) ezfio.set_qmcpack_qmc_coef(qph5['qmcpack/qmc_coef'][()].tolist()) ezfio.set_qmcpack_qmc_expo(qph5['qmcpack/qmc_expo'][()].tolist()) - + ezfio.set_qmcpack_qmc_pbc(qph5['qmcpack'].attrs['PBC']) ezfio.set_qmcpack_qmc_cart(qph5['qmcpack'].attrs['cart']) ezfio.set_qmcpack_qmc_pseudo(qph5['qmcpack'].attrs['Pseudo']) @@ -466,15 +486,15 @@ def convert_kpts(filename,qph5path,qmcpack=True): print("to create ezfio without qmcpack data, use 'qp_convert_h5_to_ezfio --noqmc'") raise - - + + ########################################## # # # MO Coef # # # ########################################## - - + + with h5py.File(qph5path,'r') as qph5: mo_coef_kpts = qph5['mo_basis/mo_coef_kpts'][()].tolist() mo_coef_cplx = qph5['mo_basis/mo_coef_complex'][()].tolist() @@ -482,63 +502,63 @@ def convert_kpts(filename,qph5path,qmcpack=True): ezfio.set_mo_basis_mo_coef_complex(mo_coef_cplx) #maybe fix qp so we don't need this? #ezfio.set_mo_basis_mo_coef([[i for i in range(mo_num)] * ao_num]) - - + + ########################################## # # # Integrals Mono # # # ########################################## - + with h5py.File(qph5path,'r') as qph5: if 'ao_one_e_ints' in qph5.keys(): kin_ao_reim=qph5['ao_one_e_ints/ao_integrals_kinetic_kpts'][()].tolist() ovlp_ao_reim=qph5['ao_one_e_ints/ao_integrals_overlap_kpts'][()].tolist() ne_ao_reim=qph5['ao_one_e_ints/ao_integrals_n_e_kpts'][()].tolist() - + ezfio.set_ao_one_e_ints_ao_integrals_kinetic_kpts(kin_ao_reim) ezfio.set_ao_one_e_ints_ao_integrals_overlap_kpts(ovlp_ao_reim) ezfio.set_ao_one_e_ints_ao_integrals_n_e_kpts(ne_ao_reim) - + 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') - - + + with h5py.File(qph5path,'r') as qph5: if 'mo_one_e_ints' in qph5.keys(): kin_mo_reim=qph5['mo_one_e_ints/mo_integrals_kinetic_kpts'][()].tolist() ovlp_mo_reim=qph5['mo_one_e_ints/mo_integrals_overlap_kpts'][()].tolist() ne_mo_reim=qph5['mo_one_e_ints/mo_integrals_n_e_kpts'][()].tolist() - + ezfio.set_mo_one_e_ints_mo_integrals_kinetic_kpts(kin_mo_reim) ezfio.set_mo_one_e_ints_mo_integrals_overlap_kpts(ovlp_mo_reim) #ezfio.set_mo_one_e_ints_mo_integrals_n_e_complex(ne_mo_reim) ezfio.set_mo_one_e_ints_mo_integrals_n_e_kpts(ne_mo_reim) - + ezfio.set_mo_one_e_ints_io_mo_integrals_kinetic('Read') ezfio.set_mo_one_e_ints_io_mo_integrals_overlap('Read') #ezfio.set_mo_one_e_ints_io_mo_integrals_n_e('Read') ezfio.set_mo_one_e_ints_io_mo_integrals_n_e('Read') - + ########################################## # # # k-points # # # ########################################## - + with h5py.File(qph5path,'r') as qph5: kconserv = qph5['nuclei/kconserv'][()].tolist() - + ezfio.set_nuclei_kconserv(kconserv) ezfio.set_nuclei_io_kconserv('Read') - + ########################################## # # # Integrals Bi # # # ########################################## - + # should this be in ao_basis? ao_two_e_ints? with h5py.File(qph5path,'r') as qph5: if 'ao_two_e_ints' in qph5.keys(): @@ -552,7 +572,7 @@ def convert_kpts(filename,qph5path,qmcpack=True): dfao_reim=qph5['ao_two_e_ints/df_ao_integrals'][()].tolist() ezfio.set_ao_two_e_ints_df_ao_integrals_complex(dfao_reim) 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'] ezfio.set_ao_two_e_ints_df_num(df_num) @@ -563,13 +583,13 @@ def convert_kpts(filename,qph5path,qmcpack=True): dfmo_reim=qph5['mo_two_e_ints/df_mo_integrals'][()].tolist() ezfio.set_mo_two_e_ints_df_mo_integrals_complex(dfmo_reim) ezfio.set_mo_two_e_ints_io_df_mo_integrals('Read') - + return def convert_cplx(filename,qph5path): ezfio.set_file(filename) ezfio.set_nuclei_is_complex(True) - + with h5py.File(qph5path,'r') as qph5: kpt_num = qph5['nuclei'].attrs['kpt_num'] nucl_num = qph5['nuclei'].attrs['nucl_num'] @@ -577,25 +597,25 @@ def convert_cplx(filename,qph5path): mo_num = qph5['mo_basis'].attrs['mo_num'] elec_alpha_num = qph5['electrons'].attrs['elec_alpha_num'] elec_beta_num = qph5['electrons'].attrs['elec_beta_num'] - + 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) - + # don't multiply nuclei by kpt_num # work in k-space, not in equivalent supercell - nucl_num_per_kpt = nucl_num + nucl_num_per_kpt = nucl_num ezfio.set_nuclei_nucl_num(nucl_num_per_kpt) - + # these are totals (kpt_num * num_per_kpt) # need to change if we want to truncate orbital space within pyscf ezfio.set_ao_basis_ao_num(ao_num) ezfio.set_mo_basis_mo_num(mo_num) ezfio.electrons_elec_alpha_num = elec_alpha_num ezfio.electrons_elec_beta_num = elec_beta_num - - - + + + ##ao_num = mo_num ##Important ! #import math @@ -605,56 +625,56 @@ def convert_cplx(filename,qph5path): # #ezfio.electrons_elec_alpha_num = int(nelec_alpha_per_kpt * n_kpts) #ezfio.electrons_elec_beta_num = int(nelec_beta_per_kpt * n_kpts) - + #ezfio.electrons_elec_alpha_num = int(math.ceil(num_elec / 2.)) #ezfio.electrons_elec_beta_num = int(math.floor(num_elec / 2.)) - + #ezfio.set_utils_num_kpts(n_kpts) #ezfio.set_integrals_bielec_df_num(n_aux) - + #(old)Important #ezfio.set_nuclei_nucl_num(nucl_num) #ezfio.set_nuclei_nucl_charge([0.]*nucl_num) #ezfio.set_nuclei_nucl_coord( [ [0.], [0.], [0.] ]*nucl_num ) #ezfio.set_nuclei_nucl_label( ['He'] * nucl_num ) - - + + with h5py.File(qph5path,'r') as qph5: nucl_charge=qph5['nuclei/nucl_charge'][()].tolist() nucl_coord=qph5['nuclei/nucl_coord'][()].T.tolist() nucl_label=qph5['nuclei/nucl_label'][()].tolist() nuclear_repulsion = qph5['nuclei'].attrs['nuclear_repulsion'] - + ezfio.set_nuclei_nucl_charge(nucl_charge) ezfio.set_nuclei_nucl_coord(nucl_coord) if isinstance(nucl_label[0],bytes): nucl_label = list(map(lambda x:x.decode(),nucl_label)) ezfio.set_nuclei_nucl_label(nucl_label) - + ezfio.set_nuclei_io_nuclear_repulsion('Read') ezfio.set_nuclei_nuclear_repulsion(nuclear_repulsion) - - + + ########################################## # # # Basis # # # ########################################## - + # with h5py.File(qph5path,'r') as qph5: # ezfio.set_ao_basis_ao_basis(qph5['ao_basis'].attrs['ao_basis']) # ezfio.set_ao_basis_ao_nucl(qph5['ao_basis/ao_nucl'][()].tolist()) -# -# +# +# # #Just need one (can clean this up later) # ao_prim_num_max = 5 -# +# # d = [ [0] *ao_prim_num_max]*ao_num # ezfio.set_ao_basis_ao_prim_num([ao_prim_num_max]*ao_num) # ezfio.set_ao_basis_ao_power(d) # ezfio.set_ao_basis_ao_coef(d) # ezfio.set_ao_basis_ao_expo(d) - + ########################################## # # # Basis # @@ -692,78 +712,78 @@ def convert_cplx(filename,qph5path): ezfio.set_ao_basis_ao_coef(qph5['ao_basis/ao_coef'][()].tolist()) ezfio.set_ao_basis_ao_expo(qph5['ao_basis/ao_expo'][()].tolist()) - - - + + + ########################################## # # # MO Coef # # # ########################################## - - + + with h5py.File(qph5path,'r') as qph5: mo_coef_reim = qph5['mo_basis/mo_coef_complex'][()].tolist() ezfio.set_mo_basis_mo_coef_complex(mo_coef_reim) #maybe fix qp so we don't need this? #ezfio.set_mo_basis_mo_coef([[i for i in range(mo_num)] * ao_num]) - - + + ########################################## # # # Integrals Mono # # # ########################################## - + with h5py.File(qph5path,'r') as qph5: if 'ao_one_e_ints' in qph5.keys(): kin_ao_reim=qph5['ao_one_e_ints/ao_integrals_kinetic'][()].tolist() ovlp_ao_reim=qph5['ao_one_e_ints/ao_integrals_overlap'][()].tolist() ne_ao_reim=qph5['ao_one_e_ints/ao_integrals_n_e'][()].tolist() - + ezfio.set_ao_one_e_ints_ao_integrals_kinetic_complex(kin_ao_reim) ezfio.set_ao_one_e_ints_ao_integrals_overlap_complex(ovlp_ao_reim) ezfio.set_ao_one_e_ints_ao_integrals_n_e_complex(ne_ao_reim) - + 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') - - + + with h5py.File(qph5path,'r') as qph5: if 'mo_one_e_ints' in qph5.keys(): kin_mo_reim=qph5['mo_one_e_ints/mo_integrals_kinetic'][()].tolist() #ovlp_mo_reim=qph5['mo_one_e_ints/mo_integrals_overlap'][()].tolist() ne_mo_reim=qph5['mo_one_e_ints/mo_integrals_n_e'][()].tolist() - + ezfio.set_mo_one_e_ints_mo_integrals_kinetic_complex(kin_mo_reim) #ezfio.set_mo_one_e_ints_mo_integrals_overlap_complex(ovlp_mo_reim) #ezfio.set_mo_one_e_ints_mo_integrals_n_e_complex(ne_mo_reim) ezfio.set_mo_one_e_ints_mo_integrals_n_e_complex(ne_mo_reim) - + ezfio.set_mo_one_e_ints_io_mo_integrals_kinetic('Read') #ezfio.set_mo_one_e_ints_io_mo_integrals_overlap('Read') #ezfio.set_mo_one_e_ints_io_mo_integrals_n_e('Read') ezfio.set_mo_one_e_ints_io_mo_integrals_n_e('Read') - + ########################################## # # # k-points # # # ########################################## - + with h5py.File(qph5path,'r') as qph5: kconserv = qph5['nuclei/kconserv'][()].tolist() - + ezfio.set_nuclei_kconserv(kconserv) ezfio.set_nuclei_io_kconserv('Read') - + ########################################## # # # Integrals Bi # # # ########################################## - + # should this be in ao_basis? ao_two_e_ints? with h5py.File(qph5path,'r') as qph5: if 'ao_two_e_ints' in qph5.keys(): @@ -777,7 +797,7 @@ def convert_cplx(filename,qph5path): dfao_reim=qph5['ao_two_e_ints/df_ao_integrals'][()].tolist() ezfio.set_ao_two_e_ints_df_ao_integrals_complex(dfao_reim) 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'] ezfio.set_ao_two_e_ints_df_num(df_num) @@ -788,9 +808,9 @@ def convert_cplx(filename,qph5path): dfmo_reim=qph5['mo_two_e_ints/df_mo_integrals'][()].tolist() ezfio.set_mo_two_e_ints_df_mo_integrals_complex(dfmo_reim) ezfio.set_mo_two_e_ints_io_df_mo_integrals('Read') - + return - + if __name__ == '__main__': ARGUMENTS = docopt(__doc__) @@ -827,4 +847,3 @@ if __name__ == '__main__': # #to be sure your MOs will be orthogonal, which is not the case when #the MOs are read from output files (not enough precision in output).""") - From 9d9cb32fd6a7f1cad012ce8983b53c8d29256053 Mon Sep 17 00:00:00 2001 From: amandadumi Date: Wed, 18 May 2022 14:50:03 -0400 Subject: [PATCH 14/28] minusk and qktok2 parse --- bin/qp_convert_h5_to_ezfio | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/bin/qp_convert_h5_to_ezfio b/bin/qp_convert_h5_to_ezfio index f1782867..b2c3a4a3 100755 --- a/bin/qp_convert_h5_to_ezfio +++ b/bin/qp_convert_h5_to_ezfio @@ -232,11 +232,13 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True): #TODO with h5py.File(qph5path,'r') as qph5: kconserv = qph5['nuclei/kconserv'][()].tolist() + minusk = qph5['Hamiltonian']['MinusK'][:]+1 + QKTok2 = qph5['Hamiltonian']['QPTok2'][:]+1 ezfio.set_nuclei_kconserv(kconserv) ezfio.set_nuclei_io_kconserv('Read') - # qktok2 - # minuxk + ezfio.set_nuclei_minusk(minusk) + ezfio.set_nuclei_qktok2(QKTok2) # kpt_sparse_map # unique_kpt_num # io_kpt_symm @@ -778,6 +780,7 @@ def convert_cplx(filename,qph5path): ezfio.set_nuclei_kconserv(kconserv) ezfio.set_nuclei_io_kconserv('Read') + ########################################## # # # Integrals Bi # From 99ce2d7418812f6f3a245d4f106f0e03792314d7 Mon Sep 17 00:00:00 2001 From: Shiv Upadhyay Date: Wed, 18 May 2022 15:11:12 -0400 Subject: [PATCH 15/28] Set unique kpt num based on number of Ls availible in hdf5 file --- bin/qp_convert_h5_to_ezfio | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/bin/qp_convert_h5_to_ezfio b/bin/qp_convert_h5_to_ezfio index b2c3a4a3..b43fb2c9 100755 --- a/bin/qp_convert_h5_to_ezfio +++ b/bin/qp_convert_h5_to_ezfio @@ -234,11 +234,15 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True): kconserv = qph5['nuclei/kconserv'][()].tolist() minusk = qph5['Hamiltonian']['MinusK'][:]+1 QKTok2 = qph5['Hamiltonian']['QPTok2'][:]+1 + unique_kpt_num = 1 + while f'L{i}' in qph5['Hamiltonian']['KPFactorized'].keys(): + unique_kpt_num += 1 ezfio.set_nuclei_kconserv(kconserv) ezfio.set_nuclei_io_kconserv('Read') ezfio.set_nuclei_minusk(minusk) ezfio.set_nuclei_qktok2(QKTok2) + ezfio.set_nuclei_unique_kpt_num(unique_kpt_num) # kpt_sparse_map # unique_kpt_num # io_kpt_symm From 0911a1c4def4728d97f7ee3dca3768c701ba8c93 Mon Sep 17 00:00:00 2001 From: Shiv Upadhyay Date: Wed, 18 May 2022 15:17:50 -0400 Subject: [PATCH 16/28] Remove loop --- bin/qp_convert_h5_to_ezfio | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/bin/qp_convert_h5_to_ezfio b/bin/qp_convert_h5_to_ezfio index b43fb2c9..1e15dab9 100755 --- a/bin/qp_convert_h5_to_ezfio +++ b/bin/qp_convert_h5_to_ezfio @@ -234,9 +234,7 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True): kconserv = qph5['nuclei/kconserv'][()].tolist() minusk = qph5['Hamiltonian']['MinusK'][:]+1 QKTok2 = qph5['Hamiltonian']['QPTok2'][:]+1 - unique_kpt_num = 1 - while f'L{i}' in qph5['Hamiltonian']['KPFactorized'].keys(): - unique_kpt_num += 1 + unique_kpt_num = len(qph5['Hamiltonian']['KPFactorized']) ezfio.set_nuclei_kconserv(kconserv) ezfio.set_nuclei_io_kconserv('Read') From 07c773b813e13f7d6aaf80f6750aaff732bd0ec3 Mon Sep 17 00:00:00 2001 From: amandadumi Date: Wed, 18 May 2022 16:28:45 -0400 Subject: [PATCH 17/28] add kpt_sparse_map --- bin/qp_convert_h5_to_ezfio | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/bin/qp_convert_h5_to_ezfio b/bin/qp_convert_h5_to_ezfio index 1e15dab9..4e7367e8 100755 --- a/bin/qp_convert_h5_to_ezfio +++ b/bin/qp_convert_h5_to_ezfio @@ -235,11 +235,20 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True): minusk = qph5['Hamiltonian']['MinusK'][:]+1 QKTok2 = qph5['Hamiltonian']['QPTok2'][:]+1 unique_kpt_num = len(qph5['Hamiltonian']['KPFactorized']) - + unique_k_idx = [] + for i in qph5['Hamiltonian']['KPFactorized'].keys(): + unique_k_idx.append(int(i[1:])+1) + kpt_sparse_map = np.zeros(kpt_num) + for i in range(kpt_num): + if i+1 is in uniq_k_idx: + kpt_sparse_map[i] = i+1 + else: + kpt_sparse_map[i] = -minusk[i] ezfio.set_nuclei_kconserv(kconserv) ezfio.set_nuclei_io_kconserv('Read') ezfio.set_nuclei_minusk(minusk) ezfio.set_nuclei_qktok2(QKTok2) + ezfio.set_nuclei_kpt_sparse_map(kpt_sparse_map) ezfio.set_nuclei_unique_kpt_num(unique_kpt_num) # kpt_sparse_map # unique_kpt_num From 272f7fdc3091ac3c615f7af7cf6b90a04f885b20 Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Wed, 29 Jun 2022 13:19:29 -0500 Subject: [PATCH 18/28] updated cd irp --- src/ao_two_e_ints/EZFIO.cfg | 5 ++--- src/ao_two_e_ints/cd_ao_ints.irp.f | 14 +++++++++++--- src/mo_two_e_ints/cd_mo_ints.irp.f | 6 +++--- 3 files changed, 16 insertions(+), 9 deletions(-) diff --git a/src/ao_two_e_ints/EZFIO.cfg b/src/ao_two_e_ints/EZFIO.cfg index 0dcb769e..aab41e1f 100644 --- a/src/ao_two_e_ints/EZFIO.cfg +++ b/src/ao_two_e_ints/EZFIO.cfg @@ -39,13 +39,12 @@ interface: ezfio type: integer doc: number of cholesky vecs for each kpt size: (nuclei.unique_kpt_num) -interface: ezfio +interface: ezfio, provider [chol_num_max] type: integer doc: max number of cholesky vecs -default: =maxval(ao_two_e_ints.chol_num) -interface: ezfio +interface: ezfio, provider [io_chol_ao_integrals] type: Disk_access diff --git a/src/ao_two_e_ints/cd_ao_ints.irp.f b/src/ao_two_e_ints/cd_ao_ints.irp.f index 866117bf..318cf6be 100644 --- a/src/ao_two_e_ints/cd_ao_ints.irp.f +++ b/src/ao_two_e_ints/cd_ao_ints.irp.f @@ -1,4 +1,12 @@ -BEGIN_PROVIDER [complex*16, chol_ao_integrals_complex, (ao_num_per_kpt,ao_num_per_kpt,chol_num_max,kpt_num,chol_unique_kpt_num)] +!BEGIN_PROVIDER [ integer, chol_num_max ] +! implicit none +! BEGIN_DOC +! ! Max number of cholesky vectors. +! END_DOC +! chol_num_max = maxval(chol_num) +!END_PROVIDER + +BEGIN_PROVIDER [complex*16, chol_ao_integrals_complex, (ao_num_per_kpt,ao_num_per_kpt,chol_num_max,kpt_num,unique_kpt_num)] implicit none BEGIN_DOC ! CD AO integrals @@ -37,7 +45,7 @@ subroutine ao_map_fill_from_chol integer :: ki,kk,kj,kl integer :: ii,ik,ij,il integer :: kikk2,kjkl2,jl2,ik2 - integer :: i_ao,j_ao,i_cd + integer :: i_ao,j_ao,i_cd,kq complex*16,allocatable :: ints_ik(:,:,:), ints_jl(:,:,:), ints_ikjl(:,:,:,:) @@ -101,7 +109,7 @@ subroutine ao_map_fill_from_chol !$OMP idx_tmp, tmp_re, tmp_im, integral,sign,use_map1) & !$OMP DEFAULT(NONE) & !$OMP SHARED(size_buffer, kpt_num, ao_num_per_kpt, ao_num_kpt_2, & - !$OMP chol_num_max, chol_num, chol_unique_kpt_num, kpt_sparse_map, qktok2, minusk, & + !$OMP chol_num_max, chol_num, unique_kpt_num, kpt_sparse_map, qktok2, minusk, & !$OMP kl,kj,kjkl2,ints_jl, & !$OMP kconserv, chol_ao_integrals_complex, ao_integrals_threshold, ao_integrals_map, ao_integrals_map_2) diff --git a/src/mo_two_e_ints/cd_mo_ints.irp.f b/src/mo_two_e_ints/cd_mo_ints.irp.f index f270397d..19f6ccb3 100644 --- a/src/mo_two_e_ints/cd_mo_ints.irp.f +++ b/src/mo_two_e_ints/cd_mo_ints.irp.f @@ -1,4 +1,4 @@ -BEGIN_PROVIDER [complex*16, chol_mo_integrals_complex, (mo_num_per_kpt,mo_num_per_kpt,chol_num_max,kpt_num,chol_unique_kpt_num)] +BEGIN_PROVIDER [complex*16, chol_mo_integrals_complex, (mo_num_per_kpt,mo_num_per_kpt,chol_num_max,kpt_num,unique_kpt_num)] implicit none BEGIN_DOC ! CD MO integrals @@ -10,7 +10,7 @@ BEGIN_PROVIDER [complex*16, chol_mo_integrals_complex, (mo_num_per_kpt,mo_num_pe print *, 'CD MO integrals read from disk' else call chol_mo_from_chol_ao(chol_mo_integrals_complex,chol_ao_integrals_complex,mo_num_per_kpt,ao_num_per_kpt, & - chol_num_max,kpt_num,chol_unique_kpt_num) + chol_num_max,kpt_num,unique_kpt_num) endif if (write_chol_mo_integrals) then @@ -126,7 +126,7 @@ subroutine mo_map_fill_from_chol_dot !$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, chol_num, mo_num_per_kpt, mo_num_kpt_2, & + !$OMP SHARED(size_buffer, kpt_num, mo_num_per_kpt, mo_num_kpt_2, & !$OMP kl,kj,kjkl2,ints_jl, & !$OMP ki,kk,kikk2,ints_ik, & !$OMP kQ, Q_idx, chol_num, & From 7b3b10e14347c3acda07abf78d661ae64006a886 Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Wed, 29 Jun 2022 13:20:22 -0500 Subject: [PATCH 19/28] working on rmg converter --- bin/qp_convert_h5_to_ezfio | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/bin/qp_convert_h5_to_ezfio b/bin/qp_convert_h5_to_ezfio index 9b6910af..04a803a2 100755 --- a/bin/qp_convert_h5_to_ezfio +++ b/bin/qp_convert_h5_to_ezfio @@ -3,13 +3,13 @@ convert hdf5 output (e.g. from PySCF) to ezfio Usage: - qp_convert_h5_to_ezfio [--noqmc] [-o EZFIO_DIR] FILE + qp_convert_h5_to_ezfio [--noqmc] [--rmg] [-o EZFIO_DIR] FILE Options: -o --output=EZFIO_DIR Produced directory by default is FILE.ezfio --noqmc don't include basis, cell, etc. for QMCPACK - -cd h5 contains cholesky decomposition informatin, these h5 result from RMG and the pyscf AFQMC converter of QMCPACK. + --rmg h5 contains cholesky decomposition informatin, these h5 result from RMG and the pyscf AFQMC converter of QMCPACK. """ from ezfio import ezfio @@ -175,6 +175,7 @@ def convert_mol(filename,qph5path): return def convert_kpts_cd(filename,qph5path,qmcpack=True): + import json ezfio.set_file(filename) ezfio.set_nuclei_is_complex(True) @@ -275,7 +276,7 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True): unique_k_idx.append(int(i[1:])+1) kpt_sparse_map = np.zeros(kpt_num) for i in range(kpt_num): - if i+1 is in uniq_k_idx: + if i+1 in uniq_k_idx: kpt_sparse_map[i] = i+1 else: kpt_sparse_map[i] = -minusk[i] @@ -295,7 +296,7 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True): ########################################## # should this be in ao_basis? ao_two_e_ints? - with h5py.File(qph5path,'r') as qph5: + with h5py.File(qph5path,'r') as qph5: nchol_per_kpt = qph5['Hamiltonian']['NCholPerKP'][:] nchol_per_kpt = nchol_per_kpt[nchol_per_kpt != 0] nchol_per_kpt_max = max(nchol_per_kpt) @@ -871,6 +872,7 @@ if __name__ == '__main__': FILE = get_full_path(ARGUMENTS['FILE']) qmcpack = True + rmg = False if ARGUMENTS["--output"]: EZFIO_FILE = get_full_path(ARGUMENTS["--output"]) else: @@ -881,10 +883,13 @@ if __name__ == '__main__': if ARGUMENTS["--rmg"]: rmg = True with h5py.File(FILE,'r') as qph5: - do_kpts = ('kconserv' in qph5['nuclei'].keys()) - if (do_kpts): + try: + do_kpts = ('kconserv' in qph5['nuclei'].keys()) + except: + do_kpts = False + if (do_kpts or rmg): print("converting HDF5 to EZFIO for periodic system") - if cd: + if rmg: print("Using RMG and AFQMC h5") convert_kpts_cd(EZFIO_FILE,FILE,qmcpack) else: From 1a9a6b240d014c0cf6cd3a5295e6d01aba0ba4c4 Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Wed, 24 Aug 2022 14:19:55 -0500 Subject: [PATCH 20/28] working on pw integral converter --- bin/qp_convert_h5_to_ezfio | 135 ++++++++++++++++++++++++++++--------- 1 file changed, 103 insertions(+), 32 deletions(-) diff --git a/bin/qp_convert_h5_to_ezfio b/bin/qp_convert_h5_to_ezfio index 04a803a2..611c0066 100755 --- a/bin/qp_convert_h5_to_ezfio +++ b/bin/qp_convert_h5_to_ezfio @@ -22,12 +22,29 @@ import gzip #fname = sys.argv[1] #qph5name = sys.argv[2] +def kconserv_p_from_qkk2_mk(qkk2,mk): + nk, nk2 = qkk2.shape + assert(nk == nk2) + kcon_p = np.zeros((nk,nk,nk),dtype=int) + for i in range(nk): + for j in range(nk): + for k in range(nk): + kcon_p[i,j,k] = qkk2[mk[j],qkk2[k,i]] + return kcon_p + + def get_full_path(file_path): file_path = os.path.expanduser(file_path) file_path = os.path.expandvars(file_path) # file_path = os.path.abspath(file_path) return file_path +def make_reim_identity_kblocks(nk,nm,na=None): + if na is None: + na = nm + single_block = np.eye(nm, na, dtype=np.complex128).view(dtype=np.float64).reshape((nm, na, 2)) + kblocks = np.tile(single_block,[nk, 1, 1, 1]) + return kblocks def flatten(l): res = [] @@ -174,7 +191,7 @@ def convert_mol(filename,qph5path): return -def convert_kpts_cd(filename,qph5path,qmcpack=True): +def convert_kpts_cd(filename,qph5path,qmcpack=True,is_ao=True): import json ezfio.set_file(filename) @@ -186,20 +203,36 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True): ezfio.set_nuclei_nucl_coord( [ [0.], [0.], [0.] ]*nucl_num ) ezfio.set_nuclei_nucl_label( ['He'] * nucl_num ) with h5py.File(qph5path,'r') as qph5: - kpt_num = qph5['Hamiltonian']['KPoints'].shape[0] - elec_alpha_num = qph5['Hamiltonian']['dims'][4] - elec_beta_num = qph5['Hamiltonian']['dims'][5] - orb_num = qph5['Hamiltonian']['dims'][3] - try: - is_ao = json.loads(qph5['metadata'][()].decode("utf-8").replace("'",'"'))['ortho_ao'] - if is_ao: - ao_num = orb_num - elif is_ao ==False: - mo_num = orb_num - else: - raise ValueError('Problem with ortho_ao key in metadata') - except: - raise UnicodeDecodeError('metadata not correctly parsed from HDF5 file') + kpt_num = qph5['Hamiltonian/KPoints'][()].shape[0] + ham_dims = qph5['Hamiltonian/dims'][()] + NMOPerKP = qph5['Hamiltonian/NMOPerKP'][()] + _, _, kpt_num, orb_num, elec_alpha_num_tot, elec_beta_num_tot, _, nchol_maybe = ham_dims + + #for now, all kpts must have same number of MOs + for nmoi in NMOPerKP: + if nmoi != NMOPerKP[0]: + print("ERROR: all KPs must have same number of MOs") + raise ValueError + + #TODO: fix na, nb in rmg + assert(elec_alpha_num_tot % kpt_num == 0) + assert(elec_beta_num_tot % kpt_num == 0) + elec_alpha_num_per_kpt = elec_alpha_num_tot // kpt_num + elec_beta_num_per_kpt = elec_beta_num_tot // kpt_num + #elec_alpha_num_per_kpt = qph5['Hamiltonian']['dims'][4] + #elec_beta_num_per_kpt = qph5['Hamiltonian']['dims'][5] + #orb_num = qph5['Hamiltonian']['dims'][3] + #try: + # is_ao = json.loads(qph5['metadata'][()].decode("utf-8").replace("'",'"'))['ortho_ao'] + # if is_ao: + # ao_num = orb_num + # elif is_ao ==False: + # mo_num = orb_num + # else: + # raise ValueError('Problem with ortho_ao key in metadata') + #except: + # raise UnicodeDecodeError('metadata not correctly parsed from HDF5 file') + ezfio.set_nuclei_kpt_num(kpt_num) kpt_pair_num = (kpt_num*kpt_num + kpt_num)//2 @@ -210,12 +243,20 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True): ezfio.set_nuclei_nucl_num(nucl_num_per_kpt) # these are totals (kpt_num * num_per_kpt) # need to change if we want to truncate orbital space within pyscf - ezfio.set_ao_basis_ao_num(ao_num) - ezfio.set_mo_basis_mo_num(mo_num) - ezfio.set_ao_basis_ao_num_per_kpt(ao_num//kpt_num) - ezfio.set_mo_basis_mo_num_per_kpt(mo_num//kpt_num) - ezfio.electrons_elec_alpha_num = elec_alpha_num - ezfio.electrons_elec_beta_num = elec_beta_num + #if is_ao: + # ao_num = orb_num*kpt_num + #TODO: fix this? + ao_num_tot = orb_num + ao_num_per_kpt = NMOPerKP[0] + mo_num_tot = orb_num + mo_num_per_kpt = NMOPerKP[0] + #mo_num_per_kpt = ao_num_per_kpt + ezfio.set_ao_basis_ao_num(ao_num_per_kpt * kpt_num) + ezfio.set_mo_basis_mo_num(mo_num_per_kpt * kpt_num) + ezfio.set_ao_basis_ao_num_per_kpt(ao_num_per_kpt) + ezfio.set_mo_basis_mo_num_per_kpt(mo_num_per_kpt) + ezfio.electrons_elec_alpha_num = elec_alpha_num_per_kpt * kpt_num + ezfio.electrons_elec_beta_num = elec_beta_num_per_kpt * kpt_num ########################################## # # # Basis # @@ -228,13 +269,38 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True): # # ########################################## #TODO + #coef_per_kpt = np.eye(mo_num_per_kpt, ao_num_per_kpt, dtype=np.complex128).view(dtype=np.float64).reshape((mo_num_per_kpt, ao_num_per_kpt, 2)) + + #mo_coef_kpts = np.tile(coef_per_kpt,[kpt_num, 1, 1, 1]) + #qph5.create_dataset('mo_basis/mo_coef_kpts',data=make_reim_identity_kblocks(kpt_num, mo_num_per_kpt, ao_num_per_kpt)) + ezfio.set_mo_basis_mo_coef_kpts(make_reim_identity_kblocks(kpt_num, mo_num_per_kpt, ao_num_per_kpt)) + ########################################## # # # Integrals Mono # # # ########################################## + + with h5py.File(qph5path,'r') as qph5: + # we don't have separate kinetic, nuc-elec, pseudo 1e ints, so just combine in nuc-elec and set rest to zero + mono_ints_tot = np.zeros((kpt_num,ao_num_per_kpt,ao_num_per_kpt,2),dtype=np.float64) + for i in range(kpt_num): + mono_ints_tot[i] = qph5[f'Hamiltonian/H1_kp{i}'][()] + ovlp_ao_reim = make_reim_identity_kblocks(kpt_num,ao_num_per_kpt,ao_num_per_kpt) + kin_ao_reim = np.zeros((kpt_num,ao_num_per_kpt,ao_num_per_kpt,2),dtype=np.float64) + ne_ao_reim = mono_ints_tot + + ezfio.set_ao_one_e_ints_ao_integrals_kinetic_kpts(kin_ao_reim) + ezfio.set_ao_one_e_ints_ao_integrals_overlap_kpts(ovlp_ao_reim) + ezfio.set_ao_one_e_ints_ao_integrals_n_e_kpts(ne_ao_reim) + + 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') + + """ - with h5py.File(qph5path,'r') as qph5: + with h5py.File(qph5path,'r') as qph5: if is_ao: kin_ao_reim= ovlp_ao_reim= @@ -267,16 +333,19 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True): ########################################## #TODO with h5py.File(qph5path,'r') as qph5: - kconserv = qph5['nuclei/kconserv'][()].tolist() - minusk = qph5['Hamiltonian']['MinusK'][:]+1 - QKTok2 = qph5['Hamiltonian']['QPTok2'][:]+1 + #kconserv = qph5['nuclei/kconserv'][()].tolist() + #TODO: change this after rmg is fixed + #minusk = qph5['Hamiltonian']['MinusK'][:]+1 + QKTok2 = qph5['Hamiltonian']['QKTok2'][:]+1 + minusk = QKTok2[:,0] + kconserv = kconserv_p_from_qkk2_mk(QKTok2-1,minusk-1)+1 unique_kpt_num = len(qph5['Hamiltonian']['KPFactorized']) unique_k_idx = [] for i in qph5['Hamiltonian']['KPFactorized'].keys(): unique_k_idx.append(int(i[1:])+1) kpt_sparse_map = np.zeros(kpt_num) for i in range(kpt_num): - if i+1 in uniq_k_idx: + if i+1 in unique_k_idx: kpt_sparse_map[i] = i+1 else: kpt_sparse_map[i] = -minusk[i] @@ -298,16 +367,18 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True): # should this be in ao_basis? ao_two_e_ints? with h5py.File(qph5path,'r') as qph5: nchol_per_kpt = qph5['Hamiltonian']['NCholPerKP'][:] + print(nchol_per_kpt) nchol_per_kpt = nchol_per_kpt[nchol_per_kpt != 0] + print(nchol_per_kpt) nchol_per_kpt_max = max(nchol_per_kpt) - ezfio.set_chol_num(nchol_per_kpt) - ezfio.set_chol_num_max(nchol_per_kpt_max) + ezfio.set_ao_two_e_ints_chol_num(nchol_per_kpt) + ezfio.set_ao_two_e_ints_chol_num_max(nchol_per_kpt_max) if is_ao: - ao_num_per_kpt = ao_num//kpt_num - ezfio.set_io_chol_ao_integrals('Read') + #ao_num_per_kpt = ao_num//kpt_num + ezfio.set_ao_two_e_ints_io_chol_ao_integrals('Read') #ao_chol_two_e_ints = np.zeros((2, ao_num_per_kpt, ao_num_per_kpt, nchol_per_kpt_max, kpt_num, len(nchol_per_kpt))) L_list = [] - for i in len(nchol_per_kpt): + for i in range(len(nchol_per_kpt)): L = qph5['Hamiltonian']['KPFactorized'][f'L{i}'][:] L.reshape(kpt_num, ao_num_per_kpt, ao_num_per_kpt, nchol_per_kpt[i], 2) L = np.einsum("ijklm->ilkjm", A, B) @@ -344,7 +415,7 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True): ezfio.set_mo_two_e_ints_df_mo_integrals_complex(dfmo_reim) ezfio.set_mo_two_e_ints_io_df_mo_integrals('Read') """ - mo_num_per_kpt = ao_num//kpt_num + #mo_num_per_kpt = ao_num//kpt_num ezfio.set_io_chol_mo_integrals('Read') #ao_chol_two_e_ints = np.zeros((2, ao_num_per_kpt, ao_num_per_kpt, nchol_per_kpt_max, kpt_num, len(nchol_per_kpt))) L_list = [] From 0335f08082c2da834ad4ea3eb203ef9e466b3a40 Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Wed, 24 Aug 2022 16:39:42 -0500 Subject: [PATCH 21/28] cholesky converter fixes --- bin/qp_convert_h5_to_ezfio | 98 +++++++++++++++++++------ src/ao_two_e_ints/two_e_integrals.irp.f | 5 ++ 2 files changed, 79 insertions(+), 24 deletions(-) diff --git a/bin/qp_convert_h5_to_ezfio b/bin/qp_convert_h5_to_ezfio index 611c0066..5fb5554b 100755 --- a/bin/qp_convert_h5_to_ezfio +++ b/bin/qp_convert_h5_to_ezfio @@ -46,6 +46,13 @@ def make_reim_identity_kblocks(nk,nm,na=None): kblocks = np.tile(single_block,[nk, 1, 1, 1]) return kblocks +def make_reim_identity_block_diag(nk,nm,na=None): + from scipy.linalg import block_diag + kblocks = make_reim_identity_kblocks(nk,nm,na).view(dtype=np.complex128).squeeze() + kblockdiag = block_diag(*kblocks).view(dtype=np.float64).reshape((nk*nm,nk*na,2)) + print(f'kblockdiag.shape = {kblockdiag.shape}') + return kblockdiag + def flatten(l): res = [] for i in l: @@ -193,15 +200,16 @@ def convert_mol(filename,qph5path): def convert_kpts_cd(filename,qph5path,qmcpack=True,is_ao=True): import json + from scipy.linalg import block_diag ezfio.set_file(filename) ezfio.set_nuclei_is_complex(True) # Dummy atom since AFQMC h5 has no atom information - nucl_num = 1 - ezfio.set_nuclei_nucl_num(nucl_num) - ezfio.set_nuclei_nucl_charge([0.]*nucl_num) - ezfio.set_nuclei_nucl_coord( [ [0.], [0.], [0.] ]*nucl_num ) - ezfio.set_nuclei_nucl_label( ['He'] * nucl_num ) + #nucl_num = 1 + #ezfio.set_nuclei_nucl_num(nucl_num) + #ezfio.set_nuclei_nucl_charge([0.]*nucl_num) + #ezfio.set_nuclei_nucl_coord( [ [0.], [0.], [0.] ]*nucl_num ) + #ezfio.set_nuclei_nucl_label( ['He'] * nucl_num ) with h5py.File(qph5path,'r') as qph5: kpt_num = qph5['Hamiltonian/KPoints'][()].shape[0] ham_dims = qph5['Hamiltonian/dims'][()] @@ -239,8 +247,8 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True,is_ao=True): ezfio.set_nuclei_kpt_pair_num(kpt_pair_num) # don't multiply nuclei by kpt_num # work in k-space, not in equivalent supercell - nucl_num_per_kpt = nucl_num - ezfio.set_nuclei_nucl_num(nucl_num_per_kpt) + #nucl_num_per_kpt = nucl_num + #ezfio.set_nuclei_nucl_num(nucl_num_per_kpt) # these are totals (kpt_num * num_per_kpt) # need to change if we want to truncate orbital space within pyscf #if is_ao: @@ -263,6 +271,29 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True,is_ao=True): # # ########################################## #TODO + nucl_num = 1 + ezfio.set_nuclei_nucl_num(nucl_num) + ezfio.set_nuclei_nucl_charge([0.]*nucl_num) + ezfio.set_nuclei_nucl_coord( [ [0.], [0.], [0.] ]*nucl_num ) + ezfio.set_nuclei_nucl_label( ['He'] * nucl_num ) + nucl_num_per_kpt = nucl_num + ezfio.set_nuclei_nucl_num(nucl_num_per_kpt) + ezfio.set_nuclei_io_kpt_symm('Read') + ezfio.set_ao_basis_ao_basis("dummy basis") + + #nucleus on which each AO is centered + ao_nucl = [1 for i in range(ao_num_per_kpt)]*kpt_num + ezfio.set_ao_basis_ao_nucl(ao_nucl) + + + #Just need one (can clean this up later) + ao_prim_num_max = 5 + + d = [ [0] *ao_prim_num_max]*ao_num_tot + ezfio.set_ao_basis_ao_prim_num([ao_prim_num_max]*ao_num_tot) + ezfio.set_ao_basis_ao_power(d) + ezfio.set_ao_basis_ao_coef(d) + ezfio.set_ao_basis_ao_expo(d) ########################################## # # # MOCoeff # @@ -274,6 +305,7 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True,is_ao=True): #mo_coef_kpts = np.tile(coef_per_kpt,[kpt_num, 1, 1, 1]) #qph5.create_dataset('mo_basis/mo_coef_kpts',data=make_reim_identity_kblocks(kpt_num, mo_num_per_kpt, ao_num_per_kpt)) ezfio.set_mo_basis_mo_coef_kpts(make_reim_identity_kblocks(kpt_num, mo_num_per_kpt, ao_num_per_kpt)) + ezfio.set_mo_basis_mo_coef_complex(make_reim_identity_block_diag(kpt_num, mo_num_per_kpt, ao_num_per_kpt)) ########################################## # # @@ -334,21 +366,25 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True,is_ao=True): #TODO with h5py.File(qph5path,'r') as qph5: #kconserv = qph5['nuclei/kconserv'][()].tolist() - #TODO: change this after rmg is fixed - #minusk = qph5['Hamiltonian']['MinusK'][:]+1 + minusk = qph5['Hamiltonian']['MinusK'][:]+1 QKTok2 = qph5['Hamiltonian']['QKTok2'][:]+1 - minusk = QKTok2[:,0] + #TODO: change this after rmg is fixed + #minusk = QKTok2[:,0] kconserv = kconserv_p_from_qkk2_mk(QKTok2-1,minusk-1)+1 unique_kpt_num = len(qph5['Hamiltonian']['KPFactorized']) unique_k_idx = [] for i in qph5['Hamiltonian']['KPFactorized'].keys(): unique_k_idx.append(int(i[1:])+1) - kpt_sparse_map = np.zeros(kpt_num) + unique_k_idx.sort() + kpt_sparse_map = np.zeros(kpt_num,dtype=int) + isparse=0 + #TODO: make robust: this assumes that for each pair, the one with data has a lower index for i in range(kpt_num): if i+1 in unique_k_idx: - kpt_sparse_map[i] = i+1 + kpt_sparse_map[i] = isparse+1 + isparse += 1 else: - kpt_sparse_map[i] = -minusk[i] + kpt_sparse_map[i] = -kpt_sparse_map[minusk[i]-1] ezfio.set_nuclei_kconserv(kconserv) ezfio.set_nuclei_io_kconserv('Read') ezfio.set_nuclei_minusk(minusk) @@ -366,10 +402,18 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True,is_ao=True): # should this be in ao_basis? ao_two_e_ints? with h5py.File(qph5path,'r') as qph5: - nchol_per_kpt = qph5['Hamiltonian']['NCholPerKP'][:] - print(nchol_per_kpt) - nchol_per_kpt = nchol_per_kpt[nchol_per_kpt != 0] + nchol_per_kpt_all = qph5['Hamiltonian']['NCholPerKP'][:] + print(nchol_per_kpt_all) + #nchol_per_kpt = nchol_per_kpt_all[nchol_per_kpt_all != 0] + nchol_per_kpt = nchol_per_kpt_all[np.array(unique_k_idx,dtype=int)-1] print(nchol_per_kpt) + print(unique_k_idx) + #for i in range(kpt_num): + # if i+1 in unique_k_idx: + # print('* ',i,nchol_per_kpt_all[i]) + # else: + # print(' ',i,nchol_per_kpt_all[i]) + nchol_per_kpt_max = max(nchol_per_kpt) ezfio.set_ao_two_e_ints_chol_num(nchol_per_kpt) ezfio.set_ao_two_e_ints_chol_num_max(nchol_per_kpt_max) @@ -378,11 +422,16 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True,is_ao=True): ezfio.set_ao_two_e_ints_io_chol_ao_integrals('Read') #ao_chol_two_e_ints = np.zeros((2, ao_num_per_kpt, ao_num_per_kpt, nchol_per_kpt_max, kpt_num, len(nchol_per_kpt))) L_list = [] - for i in range(len(nchol_per_kpt)): - L = qph5['Hamiltonian']['KPFactorized'][f'L{i}'][:] - L.reshape(kpt_num, ao_num_per_kpt, ao_num_per_kpt, nchol_per_kpt[i], 2) - L = np.einsum("ijklm->ilkjm", A, B) - L_list.append(L) + L_all = np.zeros((unique_kpt_num, kpt_num, ao_num_per_kpt, ao_num_per_kpt, nchol_per_kpt_max,2),dtype=np.float64) + print(kpt_sparse_map) + print(np.array(unique_k_idx)-1) + for i in range(unique_kpt_num): + ki = unique_k_idx[i]-1 + print(i, ki) + L_i = qph5[f'Hamiltonian/KPFactorized/L{ki}'][()].reshape((kpt_num, ao_num_per_kpt, ao_num_per_kpt, nchol_per_kpt[kpt_sparse_map[ki]-1], 2)) + #L.reshape(kpt_num, ao_num_per_kpt, ao_num_per_kpt, nchol_per_kpt[i], 2) + #L = np.einsum("ijklm->ilkjm", A, B) + L_all[i,:,:,:,:nchol_per_kpt[kpt_sparse_map[ki]-1],:] = L_i #(6, 5184, 2) """ @@ -393,9 +442,10 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True,is_ao=True): for kpt_idx in range(kpt_num): ao_chol_two_e_ints[cmplx][ao_idx_i][ao_idx_j][chol_idx][kpt_idx][i] = L[kpt_idx][ao_idx_i][ao_idx_j][chol_idx][cmplx] """ - ao_chol_two_e_ints = np.vstack(L_list) - ao_chol_two_e_ints = ao_chol_two_e_ints.transpose() - ezfio.set_chol_ao_integrals_complex(ao_chol_two_e_ints) + #ao_chol_two_e_ints = np.vstack(L_list) + #ao_chol_two_e_ints = ao_chol_two_e_ints.transpose() + #TODO: check dims/reshape/transpose + ezfio.set_ao_two_e_ints_chol_ao_integrals_complex(L_all) #(2,ao_basis.ao_num_per_kpt,ao_basis.ao_num_per_kpt,ao_two_e_ints.chol_num_max,nuclei.kpt_num,nuclei.unique_kpt_num#) """ diff --git a/src/ao_two_e_ints/two_e_integrals.irp.f b/src/ao_two_e_ints/two_e_integrals.irp.f index c2a48a2f..662492b8 100644 --- a/src/ao_two_e_ints/two_e_integrals.irp.f +++ b/src/ao_two_e_ints/two_e_integrals.irp.f @@ -358,6 +358,11 @@ BEGIN_PROVIDER [ logical, ao_two_e_integrals_in_map ] print*, 'AO integrals provided from 3-index ao ints (periodic)' ao_two_e_integrals_in_map = .True. return + else if (read_chol_ao_integrals) then + call ao_map_fill_from_chol + print*, 'AO integrals provided from 3-index Cholesky ao ints (periodic)' + ao_two_e_integrals_in_map = .True. + return else print*,'calculation of periodic AOs not implemented' stop -1 From 2371798bdc8b1157439e9c69a9a269a583c8f0c9 Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Thu, 15 Sep 2022 15:34:26 -0500 Subject: [PATCH 22/28] debugging cd block gemm --- src/ao_two_e_ints/cd_ao_ints.irp.f | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/ao_two_e_ints/cd_ao_ints.irp.f b/src/ao_two_e_ints/cd_ao_ints.irp.f index 318cf6be..19305d42 100644 --- a/src/ao_two_e_ints/cd_ao_ints.irp.f +++ b/src/ao_two_e_ints/cd_ao_ints.irp.f @@ -102,7 +102,7 @@ subroutine ao_map_fill_from_chol enddo endif - !$OMP PARALLEL PRIVATE(i,k,j,l,ki,kk,ii,ik,ij,il,kikk2,jl2,ik2,kQ, & + !$OMP PARALLEL PRIVATE(i,k,j,l,ki,kk,ii,ik,ij,il,kikk2,jl2,ik2, & !$OMP ints_ik, ints_ikjl, i_ao, j_ao, i_cd, & !$OMP n_integrals_1, buffer_i_1, buffer_values_1, & !$OMP n_integrals_2, buffer_i_2, buffer_values_2, & @@ -110,7 +110,7 @@ subroutine ao_map_fill_from_chol !$OMP DEFAULT(NONE) & !$OMP SHARED(size_buffer, kpt_num, ao_num_per_kpt, ao_num_kpt_2, & !$OMP chol_num_max, chol_num, unique_kpt_num, kpt_sparse_map, qktok2, minusk, & - !$OMP kl,kj,kjkl2,ints_jl, & + !$OMP kl,kj,kjkl2,ints_jl,kQ, & !$OMP kconserv, chol_ao_integrals_complex, ao_integrals_threshold, ao_integrals_map, ao_integrals_map_2) allocate( & @@ -124,6 +124,8 @@ subroutine ao_map_fill_from_chol !$OMP DO SCHEDULE(guided) do kk=1,kl + !print*,'debug' + !print*,kQ,kl,kj,kk ki = qktok2(minusk(kk),kQ) assert(ki == kconserv(kl,kk,kj)) if (ki>kl) cycle From fe9ddc4d988acec2d344495344b288b22e355c59 Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Wed, 21 Sep 2022 15:47:51 -0500 Subject: [PATCH 23/28] debugging cd ints --- bin/qp_convert_h5_to_ezfio | 360 ++++++++++++++++++++++++++- src/ao_two_e_ints/cd_ao_ints.irp.f | 2 + src/mo_two_e_ints/cd_mo_ints.irp.f | 10 +- src/utils_complex/dump_cd_ints.irp.f | 29 +++ src/utils_complex/dump_cd_ksym.irp.f | 47 ++++ src/utils_complex/test_cd_ksym.irp.f | 234 +++++++++++++++++ 6 files changed, 675 insertions(+), 7 deletions(-) create mode 100644 src/utils_complex/dump_cd_ints.irp.f create mode 100644 src/utils_complex/dump_cd_ksym.irp.f create mode 100644 src/utils_complex/test_cd_ksym.irp.f diff --git a/bin/qp_convert_h5_to_ezfio b/bin/qp_convert_h5_to_ezfio index 5fb5554b..5c8f8ad3 100755 --- a/bin/qp_convert_h5_to_ezfio +++ b/bin/qp_convert_h5_to_ezfio @@ -21,7 +21,93 @@ from docopt import docopt import gzip #fname = sys.argv[1] #qph5name = sys.argv[2] +def idx2_tri(i,j): + """ + for 0-indexed counting + """ + p = max(i,j) + q = min(i,j) + return q + (p*(p+1))//2 +def idx2_tri_1(i,j): + """ + for 1-indexed counting + """ + p = max(i,j) + q = min(i,j) + return q + (p*(p-1))//2 +def idx4_cplx_1(i,j,k,l): + """ + original function from qp2 (fortran counting) + """ + p = idx2_tri_1(i,k) + q = idx2_tri_1(j,l) + i1 = idx2_tri_1(p,q) + return (i1,p,q) + +def ao_idx_map_sign(i,j,k,l): + """ + qp2 indexing + """ + idx,ik,jl = idx4_cplx_1(i,j,k,l) + ij = idx2_tri_1(i,j) + kl = idx2_tri_1(k,l) + idx = 2*idx - 1 + + if (ij==kl): + sign = 0.0 + use_map1 = False + else: + if ik==jl: + if i25s} {fortformat(vi):>25s}') + + if dump_fci: + Wfull = np.zeros((ao_num_tot, ao_num_tot, ao_num_tot, ao_num_tot), dtype=np.complex128) + for Qi in range(kpt_num): + Qloc = abs(kpt_sparse_map[Qi])-1 + Qneg = (kpt_sparse_map[Qi] < 0) + LQ00 = L_all[Qloc] + #LQ0a = LQ00.view(dtype=np.complex128) + #print(f'LQ0a.shape {LQ0a.shape}') + #LQ0a1 = LQ0a.reshape((kpt_num,ao_num_per_kpt,ao_num_per_kpt,-1)) + #print(f'LQ0a1.shape {LQ0a1.shape}') + LQ0 = LQ00[:,:,:,:,0] + 1j*LQ00[:,:,:,:,1] + #print(f'LQ0.shape {LQ0.shape}') + #print(f'abdiff {np.abs(LQ0a1 - LQ0).max()}') + + + + for kp in range(kpt_num): + kr = QKTok2[Qi,kp]-1 + for ks in range(kpt_num): + kq = QKTok2[Qi,ks]-1 + if Qneg: + A = LQ0[kr].transpose((1,0,2)).conj() + B = LQ0[kq] + W = np.einsum('prn,sqn->pqrs',A,B) + else: + A = LQ0[kp] + B = LQ0[ks].transpose((1,0,2)).conj() + W = np.einsum('rpn,qsn->pqrs',A,B) + p0 = kp*ao_num_per_kpt + r0 = kr*ao_num_per_kpt + q0 = kq*ao_num_per_kpt + s0 = ks*ao_num_per_kpt + for ip in range(ao_num_per_kpt): + for iq in range(ao_num_per_kpt): + for ir in range(ao_num_per_kpt): + for i_s in range(ao_num_per_kpt): + v = W[ip,iq,ir,i_s] + print(f'{p0+ip:5d} {q0+iq:5d} {r0+ir:5d} {s0+i_s:5d} {v.real:25.15E} {v.imag:25.15E}') + Wfull[p0:p0+ao_num_per_kpt,q0:q0+ao_num_per_kpt,r0:r0+ao_num_per_kpt,s0:s0+ao_num_per_kpt] = W.copy() + H1 = np.zeros((ao_num_tot, ao_num_tot), dtype=np.complex128) + for Qi in range(kpt_num): + hi0 = qph5[f'Hamiltonian/H1_kp{Qi}'][()] + hi = hi0[:,:,0] + 1j*hi0[:,:,1] + H1[Qi*ao_num_per_kpt:(Qi+1)*ao_num_per_kpt,Qi*ao_num_per_kpt:(Qi+1)*ao_num_per_kpt] = hi.copy() + mo_occ = ([1,]* elec_alpha_num_per_kpt + [0,] * (ao_num_per_kpt - elec_alpha_num_per_kpt)) * kpt_num + E1=0 + E2j=0 + E2k=0 + for imo, iocc in enumerate(mo_occ): + if iocc: + E1 += 2*H1[imo,imo] + for jmo, jocc in enumerate(mo_occ): + if jocc: + E2j += 2* Wfull[imo,jmo,imo,jmo] + E2k -= Wfull[imo,jmo,jmo,imo] + print(f'E1 = {E1:25.15E}') + print(f'E2j = {E2j:25.15E}') + print(f'E2k = {E2k:25.15E}') + print(f'E2 = {E2j+E2k:25.15E}') + + + if dump_fci2: + ao_map1 = {} + ao_map2 = {} + ao_map1_idx = {} + ao_map2_idx = {} + Wdict = {} + for kQ in range(kpt_num): + for kl in range(kpt_num): + kj = QKTok2[kQ,kl]-1 + if (kj>kl): + continue + kjkl2 = idx2_tri(kj,kl) + Qneg = (kpt_sparse_map[kQ] < 0) + Qloc = abs(kpt_sparse_map[kQ]) - 1 + if not Qneg: + ints_jl0 = L_all[Qloc,kl,:,:,:,:] + ints_jl = ints_jl0[:,:,:,0]+1j*ints_jl0[:,:,:,1] + else: + ints_jl0 = L_all[Qloc,kj,:,:,:,:] + ints_jl1 = ints_jl0[:,:,:,0]+1j*ints_jl0[:,:,:,1] + ints_jl = ints_jl1.transpose((1,0,2)).conj() + for kk in range(kl+1): + ki = QKTok2[minusk[kk]-1,kQ]-1 #TODO: check + if ki != kconserv[kl,kk,kj]-1: + print(ki,kconserv[kl,kk,kj],kl,kk,kj) + assert( ki == kconserv[kl,kk,kj]-1) #TODO: check + if (ki > kl): + continue + kikk2 = idx2_tri(ki,kk) + if not Qneg: + ints_ik0 = L_all[Qloc,ki,:,:,:] + ints_ik = ints_ik0[:,:,:,0] + 1j*ints_ik0[:,:,:,1] + else: + ints_ik0 = L_all[Qloc,kk,:,:,:] + ints_ik1 = ints_ik0[:,:,:,0]+1j*ints_ik0[:,:,:,1] + ints_ik = ints_ik1.transpose((1,0,2)).conj() + ints_jl_flat = ints_jl.reshape((ao_num_per_kpt**2, -1)) + ints_ik_flat = ints_ik.reshape((ao_num_per_kpt**2, -1)) + + ints_ikjl = np.einsum('an,bn->ab',ints_ik_flat,ints_jl_flat).reshape((ao_num_per_kpt,)*4) + + for il in range(ao_num_per_kpt): + l = il + kl*ao_num_per_kpt + for ij in range(ao_num_per_kpt): + j = ij + kj*ao_num_per_kpt + if (j>l): + break + jl2 = idx2_tri(j,l) + for ik in range(ao_num_per_kpt): + k = ik + kk*ao_num_per_kpt + if (k>l): + break + for ii in range(ao_num_per_kpt): + i = ii + ki*ao_num_per_kpt + if ((j==l) and (i>k)): + break + ik2 = idx2_tri(i,k) + if (ik2 > jl2): + break + integral = ints_ikjl[ii,ik,ij,il] + if abs(integral) < 1E-15: + continue + idx_tmp,use_map1,sign = ao_idx_map_sign(i+1,j+1,k+1,l+1) + tmp_re = integral.real + tmp_im = integral.imag + Wdict[i,j,k,l] = tmp_re + 1j*tmp_im + if use_map1: + if idx_tmp in ao_map1: + print(idx_tmp,1) + raise + ao_map1[idx_tmp] = tmp_re + ao_map1_idx[idx_tmp] = (i,j,k,l,'re') + if sign != 0.0: + ao_map1[idx_tmp+1] = tmp_im*sign + ao_map1_idx[idx_tmp+1] = (i,j,k,l,'im') + else: + if idx_tmp in ao_map2: + print(idx_tmp,2) + raise + ao_map2[idx_tmp] = tmp_re + ao_map2_idx[idx_tmp] = (i,j,k,l,'re') + if sign != 0.0: + ao_map2[idx_tmp+1] = tmp_im*sign + ao_map2_idx[idx_tmp+1] = (i,j,k,l,'im') +# for idx in ao_map1: +# i,j,k,l,ax = ao_map1_idx[idx] +# print(f'1,{idx},{i},{j},{k},{l},{ax},{ao_map1[idx]}') +# for idx in ao_map2: +# i,j,k,l,ax = ao_map2_idx[idx] +# print(f'2,{idx},{i},{j},{k},{l},{ax},{ao_map2[idx]}') + for idx in Wdict: + i,j,k,l = idx + v = Wdict[idx] + print(f'{i+1:6d} {j+1:6d} {k+1:6d} {l+1:6d} {fortformat(v.real):>25s} {fortformat(v.imag):>25s}') + + + + + + #for Qi in range(kpt_num): + # Qloc = abs(kpt_sparse_map[Qi])-1 + # Qneg = (kpt_sparse_map[Qi] < 0) + # LQ00 = L_all[Qloc] + # #LQ0a = LQ00.view(dtype=np.complex128) + # #print(f'LQ0a.shape {LQ0a.shape}') + # #LQ0a1 = LQ0a.reshape((kpt_num,ao_num_per_kpt,ao_num_per_kpt,-1)) + # #print(f'LQ0a1.shape {LQ0a1.shape}') + # LQ0 = LQ00[:,:,:,:,0] + 1j*LQ00[:,:,:,:,1] + # #print(f'LQ0.shape {LQ0.shape}') + # #print(f'abdiff {np.abs(LQ0a1 - LQ0).max()}') + + # + + # for kp in range(kpt_num): + # kr = QKTok2[Qi,kp]-1 + # for ks in range(kpt_num): + # kq = QKTok2[Qi,ks]-1 + # if Qneg: + # A = LQ0[kr].transpose((1,0,2)).conj() + # B = LQ0[kq] + # W = np.einsum('prn,sqn->pqrs',A,B) + # else: + # A = LQ0[kp] + # B = LQ0[ks].transpose((1,0,2)).conj() + # W = np.einsum('rpn,qsn->pqrs',A,B) + # p0 = kp*ao_num_per_kpt + # r0 = kr*ao_num_per_kpt + # q0 = kq*ao_num_per_kpt + # s0 = ks*ao_num_per_kpt + # for ip in range(ao_num_per_kpt): + # for iq in range(ao_num_per_kpt): + # for ir in range(ao_num_per_kpt): + # for i_s in range(ao_num_per_kpt): + # v = W[ip,iq,ir,i_s] + # print(f'{p0+ip:5d} {q0+iq:5d} {r0+ir:5d} {s0+i_s:5d} {v.real:25.15E} {v.imag:25.15E}') + # Wfull[p0:p0+ao_num_per_kpt,q0:q0+ao_num_per_kpt,r0:r0+ao_num_per_kpt,s0:s0+ao_num_per_kpt] = W.copy() + #H1 = np.zeros((ao_num_tot, ao_num_tot), dtype=np.complex128) + #for Qi in range(kpt_num): + # hi0 = qph5[f'Hamiltonian/H1_kp{Qi}'][()] + # hi = hi0[:,:,0] + 1j*hi0[:,:,1] + # H1[Qi*ao_num_per_kpt:(Qi+1)*ao_num_per_kpt,Qi*ao_num_per_kpt:(Qi+1)*ao_num_per_kpt] = hi.copy() + #mo_occ = ([1,]* elec_alpha_num_per_kpt + [0,] * (ao_num_per_kpt - elec_alpha_num_per_kpt)) * kpt_num + #E1=0 + #E2j=0 + #E2k=0 + #for imo, iocc in enumerate(mo_occ): + # if iocc: + # E1 += 2*H1[imo,imo] + # for jmo, jocc in enumerate(mo_occ): + # if jocc: + # E2j += 2* Wfull[imo,jmo,imo,jmo] + # E2k -= Wfull[imo,jmo,jmo,imo] + #print(f'E1 = {E1:25.15E}') + #print(f'E2j = {E2j:25.15E}') + #print(f'E2k = {E2k:25.15E}') + #print(f'E2 = {E2j+E2k:25.15E}') + + + #(2,ao_basis.ao_num_per_kpt,ao_basis.ao_num_per_kpt,ao_two_e_ints.chol_num_max,nuclei.kpt_num,nuclei.unique_kpt_num#) """ @@ -991,6 +1339,12 @@ def convert_cplx(filename,qph5path): if __name__ == '__main__': ARGUMENTS = docopt(__doc__) + #for i in range(1,6): + # for j in range(1,6): + # for k in range(1,6): + # for l in range(1,6): + # idx,usem1,sgn = ao_idx_map_sign(i,j,k,l) + # print(f'{i:4d} {j:4d} {k:4d} {l:4d} {str(usem1)[0]:s} {idx:6d} {sgn:5.1f}') FILE = get_full_path(ARGUMENTS['FILE']) qmcpack = True rmg = False diff --git a/src/ao_two_e_ints/cd_ao_ints.irp.f b/src/ao_two_e_ints/cd_ao_ints.irp.f index 19305d42..e1888bf9 100644 --- a/src/ao_two_e_ints/cd_ao_ints.irp.f +++ b/src/ao_two_e_ints/cd_ao_ints.irp.f @@ -92,11 +92,13 @@ subroutine ao_map_fill_from_chol !TODO: verify the kj, kl as 4th index in expressions below if (kpt_sparse_map(kQ) > 0) then ints_jl = chol_ao_integrals_complex(:,:,:,kl,kpt_sparse_map(kQ)) + !ints_jl = dconjg(chol_ao_integrals_complex(:,:,:,kl,kpt_sparse_map(kQ))) else do i_ao=1,ao_num_per_kpt do j_ao=1,ao_num_per_kpt do i_cd=1,chol_num_max ints_jl(i_ao,j_ao,i_cd) = dconjg(chol_ao_integrals_complex(j_ao,i_ao,i_cd,kj,-kpt_sparse_map(kQ))) + !ints_jl(i_ao,j_ao,i_cd) = chol_ao_integrals_complex(j_ao,i_ao,i_cd,kj,-kpt_sparse_map(kQ)) enddo enddo enddo diff --git a/src/mo_two_e_ints/cd_mo_ints.irp.f b/src/mo_two_e_ints/cd_mo_ints.irp.f index 19f6ccb3..1afc2729 100644 --- a/src/mo_two_e_ints/cd_mo_ints.irp.f +++ b/src/mo_two_e_ints/cd_mo_ints.irp.f @@ -81,7 +81,8 @@ subroutine mo_map_fill_from_chol_dot do i_mo=1,mo_num_per_kpt do j_mo=1,mo_num_per_kpt do i_cd=1,chol_num(kQ) - ints_jl(i_cd,i_mo,j_mo) = chol_mo_integrals_complex(i_mo,j_mo,i_cd,kl,Q_idx) + !ints_jl(i_cd,i_mo,j_mo) = chol_mo_integrals_complex(i_mo,j_mo,i_cd,kl,Q_idx) + ints_jl(i_cd,i_mo,j_mo) = dconjg(chol_mo_integrals_complex(i_mo,j_mo,i_cd,kl,Q_idx)) enddo enddo enddo @@ -89,7 +90,8 @@ subroutine mo_map_fill_from_chol_dot do i_mo=1,mo_num_per_kpt do j_mo=1,mo_num_per_kpt do i_cd=1,chol_num(kQ) - ints_jl(i_cd,i_mo,j_mo) = dconjg(chol_mo_integrals_complex(j_mo,i_mo,i_cd,kj,-Q_idx)) + !ints_jl(i_cd,i_mo,j_mo) = dconjg(chol_mo_integrals_complex(j_mo,i_mo,i_cd,kj,-Q_idx)) + ints_jl(i_cd,i_mo,j_mo) = chol_mo_integrals_complex(j_mo,i_mo,i_cd,kj,-Q_idx) enddo enddo enddo @@ -268,7 +270,7 @@ subroutine chol_mo_from_chol_ao(cd_mo,cd_ao,n_mo,n_ao,n_cd,n_k,n_unique_k) complex*16,allocatable :: coef_i(:,:), coef_k(:,:), ints_ik(:,:), ints_tmp(:,:) double precision :: wall_1,wall_2,cpu_1,cpu_2 - print*,'providing 3-index MO integrals from 3-index AO integrals' + print*,'providing 3-index CD MO integrals from 3-index CD AO integrals' cd_mo = 0.d0 @@ -315,7 +317,7 @@ subroutine chol_mo_from_chol_ao(cd_mo,cd_ao,n_mo,n_ao,n_cd,n_k,n_unique_k) ) call wall_time(wall_2) call cpu_time(cpu_2) - print*,' 3-idx MO provided' + print*,' 3-idx CD MO provided' 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),')' diff --git a/src/utils_complex/dump_cd_ints.irp.f b/src/utils_complex/dump_cd_ints.irp.f new file mode 100644 index 00000000..40adb2b4 --- /dev/null +++ b/src/utils_complex/dump_cd_ints.irp.f @@ -0,0 +1,29 @@ +program dump_cd_ksym + call run +end + +subroutine run + use map_module + implicit none + + integer ::q,k,n,i,j + double precision :: vr, vi + complex*16 :: v + print*,"chol_ao_integrals_complex q,k,n,i,j" + provide chol_ao_integrals_complex + do q = 1, unique_kpt_num + do k = 1, kpt_num + do n = 1, chol_num_max + do i = 1, ao_num_per_kpt + do j = 1, ao_num_per_kpt + v = chol_ao_integrals_complex(i,j,n,k,q) + vr = dble(v) + vi = dimag(v) + print '(5(I6,X),2(E25.15,X))', q, k, n, i, j, vr, vi + enddo + enddo + enddo + enddo + enddo + +end diff --git a/src/utils_complex/dump_cd_ksym.irp.f b/src/utils_complex/dump_cd_ksym.irp.f new file mode 100644 index 00000000..d5f9cfdb --- /dev/null +++ b/src/utils_complex/dump_cd_ksym.irp.f @@ -0,0 +1,47 @@ +program dump_cd_ksym + call run +end + +subroutine run + use map_module + implicit none + + integer ::i,j,k,l + integer(key_kind) :: idx + logical :: use_map1 + double precision :: sign + do i=1,5 + do j=1,5 + do k=1,5 + do l=1,5 + call ao_two_e_integral_complex_map_idx_sign(i,j,k,l,use_map1,idx,sign) + print'(4(I4,X),(L6),(I8),(F10.1))',i,j,k,l,use_map1,idx,sign + enddo + enddo + enddo + enddo + + provide qktok2 minusk kconserv + print*,'minusk' + do i=1,kpt_num + j = minusk(i) + print'(2(I4))',i,j + enddo + print*,'qktok2' + do i=1,kpt_num + do j=1,kpt_num + k = qktok2(i,j) + print'(3(I4))',i,j,k + enddo + enddo + print*,'kconserv' + do i=1,kpt_num + do j=1,kpt_num + do k=1,kpt_num + l = kconserv(i,j,k) + print'(4(I4))',i,j,k,l + enddo + enddo + enddo + +end diff --git a/src/utils_complex/test_cd_ksym.irp.f b/src/utils_complex/test_cd_ksym.irp.f new file mode 100644 index 00000000..be1433be --- /dev/null +++ b/src/utils_complex/test_cd_ksym.irp.f @@ -0,0 +1,234 @@ +program test_cd_ksym + call run +end + +subroutine run + use map_module + implicit none + + !integer ::i,j,k,l + + provide qktok2 minusk kconserv + !print*,'minusk' + !do i=1,kpt_num + ! j = minusk(i) + ! print'(2(I4))',i,j + !enddo + !print*,'qktok2' + !do i=1,kpt_num + ! do j=1,kpt_num + ! k = qktok2(i,j) + ! print'(3(I4))',i,j,k + ! enddo + !enddo + !print*,'kconserv' + !do i=1,kpt_num + ! do j=1,kpt_num + ! do k=1,kpt_num + ! l = kconserv(i,j,k) + ! print'(4(I4))',i,j,k,l + ! enddo + ! enddo + !enddo + + integer :: i,k,j,l + integer :: ki,kk,kj,kl + integer :: ii,ik,ij,il + integer :: kikk2,kjkl2,jl2,ik2 + integer :: i_ao,j_ao,i_cd,kq + + 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 :: ao_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 + + ao_num_kpt_2 = ao_num_per_kpt * ao_num_per_kpt + + size_buffer = min(ao_num_per_kpt*ao_num_per_kpt*ao_num_per_kpt,16000000) + print*, 'Providing the ao_bielec integrals from 3-index cholesky 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(ao_num_per_kpt,ao_num_per_kpt,chol_num_max)) + + wall_0 = wall_1 + ! ki + kj == kk + kl required for to be nonzero + !TODO: change loops so that we only iterate over "correct" slices (i.e. ik block is stored directly, not as conj. transp.) + ! possible cases for (ik,jl) are (+,+), (+,-), (-,+), (-,-) + ! where + is the slice used as stored, and - is the conj. transp. of the stored data + ! (+,+) and (-,-) give the same information; we should always use (+,+) + ! (+,-) and (-,+) give the same information; we should always use (+,-) + do kQ = 1, kpt_num + do kl = 1, kpt_num + kj = qktok2(kQ,kl) + assert(kQ == qktok2(kj,kl)) + if (kj>kl) cycle + call idx2_tri_int(kj,kl,kjkl2) + !TODO: verify the kj, kl as 4th index in expressions below + !if (kpt_sparse_map(kQ) > 0) then + ! ints_jl = chol_ao_integrals_complex(:,:,:,kl,kpt_sparse_map(kQ)) + !else + ! !do i_ao=1,ao_num_per_kpt + ! ! do j_ao=1,ao_num_per_kpt + ! ! do i_cd=1,chol_num_max + ! ! ints_jl(i_ao,j_ao,i_cd) = dconjg(chol_ao_integrals_complex(j_ao,i_ao,i_cd,kj,-kpt_sparse_map(kQ))) + ! ! enddo + ! ! enddo + ! !enddo + !endif + + !allocate( & + ! ints_ik(ao_num_per_kpt,ao_num_per_kpt,chol_num_max), & + ! ints_ikjl(ao_num_per_kpt,ao_num_per_kpt,ao_num_per_kpt,ao_num_per_kpt), & + ! buffer_i_1(size_buffer), & + ! buffer_i_2(size_buffer), & + ! buffer_values_1(size_buffer), & + ! buffer_values_2(size_buffer) & + !) + + do kk=1,kl + ki = qktok2(minusk(kk),kQ) + assert(ki == kconserv(kl,kk,kj)) + if (ki>kl) cycle + ! if ((kl == kj) .and. (ki > kk)) cycle + call idx2_tri_int(ki,kk,kikk2) + print*,kQ,kl,kj,kk,ki + ! if (kikk2 > kjkl2) cycle + !TODO: check this! (ki, kk slice index and transpose/notranspose) + !if (kpt_sparse_map(kQ) > 0) then + ! ints_ik = chol_ao_integrals_complex(:,:,:,ki,kpt_sparse_map(kQ)) + !else + ! do i_ao=1,ao_num_per_kpt + ! do j_ao=1,ao_num_per_kpt + ! do i_cd=1,chol_num_max + ! ints_jl(i_ao,j_ao,i_cd) = dconjg(chol_ao_integrals_complex(j_ao,i_ao,i_cd,kk,-kpt_sparse_map(kQ))) + ! enddo + ! enddo + ! enddo + !endif + + !call zgemm('N','T', ao_num_kpt_2, ao_num_kpt_2, chol_num(kQ), & + ! (1.d0,0.d0), ints_ik, ao_num_kpt_2, & + ! ints_jl, ao_num_kpt_2, & + ! (0.d0,0.d0), ints_ikjl, ao_num_kpt_2) + + !n_integrals_1=0 + !n_integrals_2=0 + !do il=1,ao_num_per_kpt + ! l=il+(kl-1)*ao_num_per_kpt + ! do ij=1,ao_num_per_kpt + ! j=ij+(kj-1)*ao_num_per_kpt + ! if (j>l) exit + ! call idx2_tri_int(j,l,jl2) + ! do ik=1,ao_num_per_kpt + ! k=ik+(kk-1)*ao_num_per_kpt + ! if (k>l) exit + ! do ii=1,ao_num_per_kpt + ! i=ii+(ki-1)*ao_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) < ao_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 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 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 insert_into_ao_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1) + !endif + !if (n_integrals_2 > 0) then + ! call insert_into_ao_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_2) + !endif + enddo !kk + !deallocate( & + ! ints_ik, & + ! ints_ikjl, & + ! buffer_i_1, & + ! buffer_i_2, & + ! buffer_values_1, & + ! buffer_values_2 & + ! ) + enddo !kl + call wall_time(wall_2) + if (wall_2 - wall_0 > 1.d0) then + wall_0 = wall_2 + !print*, 100.*float(kQ)/float(kpt_num), '% in ', & + ! wall_2-wall_1,'s',map_mb(ao_integrals_map),'+',map_mb(ao_integrals_map_2),'MB' + endif + + enddo !kQ + !deallocate( ints_jl ) + + !call map_sort(ao_integrals_map) + !call map_unique(ao_integrals_map) + !call map_sort(ao_integrals_map_2) + !call map_unique(ao_integrals_map_2) + !call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_complex_1',ao_integrals_map) + !call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_complex_2',ao_integrals_map_2) + !call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read') + + call wall_time(wall_2) + call cpu_time(cpu_2) + + !integer*8 :: get_ao_map_size, ao_map_size + !ao_map_size = get_ao_map_size() + + print*,'AO integrals provided:' + !print*,' Size of AO map ', map_mb(ao_integrals_map),'+',map_mb(ao_integrals_map_2),'MB' + !print*,' Number of AO integrals: ', ao_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 From 807a7812764ee0fe8d6c7fc734b9923bd3b6ebf3 Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Tue, 27 Sep 2022 13:49:49 -0500 Subject: [PATCH 24/28] cd fix --- bin/qp_convert_h5_to_ezfio | 36 ++++++++++++++++++++----- src/ao_two_e_ints/cd_ao_ints.irp.f | 8 +++--- src/mo_two_e_ints/cd_mo_ints.irp.f | 3 ++- src/utils_complex/dump_ao_2e_cplx.irp.f | 4 ++- 4 files changed, 39 insertions(+), 12 deletions(-) diff --git a/bin/qp_convert_h5_to_ezfio b/bin/qp_convert_h5_to_ezfio index 5c8f8ad3..4bee742f 100755 --- a/bin/qp_convert_h5_to_ezfio +++ b/bin/qp_convert_h5_to_ezfio @@ -288,8 +288,8 @@ def convert_mol(filename,qph5path): def convert_kpts_cd(filename,qph5path,qmcpack=True,is_ao=True): import json from scipy.linalg import block_diag - dump_fci, dump_cd = (False, False) - #dump_fci, dump_cd = (True, False) + #dump_fci, dump_cd = (False, False) + dump_fci, dump_cd = (True, False) #dump_fci, dump_cd = (False, True) #dump_fci2 = True dump_fci2 = False @@ -594,14 +594,33 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True,is_ao=True): kr = QKTok2[Qi,kp]-1 for ks in range(kpt_num): kq = QKTok2[Qi,ks]-1 + # 3 + #if Qneg: + # A = LQ0[kr].transpose((1,0,2)).conj() + # B = LQ0[kq] + # W = np.einsum('prn,sqn->pqrs',A,B) + #else: + # A = LQ0[kp] + # B = LQ0[ks].transpose((1,0,2)).conj() + # W = np.einsum('rpn,qsn->pqrs',A,B) + # 4 + #if Qneg: + # A = LQ0[kr].transpose((1,0,2)).conj() + # B = LQ0[kq].transpose((1,0,2)) + # W = np.einsum('prn,sqn->pqrs',A,B) + #else: + # A = LQ0[kp] + # B = LQ0[ks].conj() + # W = np.einsum('prn,sqn->pqrs',A,B) + # 5 if Qneg: A = LQ0[kr].transpose((1,0,2)).conj() - B = LQ0[kq] - W = np.einsum('prn,sqn->pqrs',A,B) + B = LQ0[kq].transpose((1,0,2)) + W = np.einsum('prn,qsn->pqrs',A,B) else: A = LQ0[kp] - B = LQ0[ks].transpose((1,0,2)).conj() - W = np.einsum('rpn,qsn->pqrs',A,B) + B = LQ0[ks].conj() + W = np.einsum('prn,qsn->pqrs',A,B) p0 = kp*ao_num_per_kpt r0 = kr*ao_num_per_kpt q0 = kq*ao_num_per_kpt @@ -622,6 +641,11 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True,is_ao=True): E1=0 E2j=0 E2k=0 + print("Jij Kij") + for i in range(ao_num_tot): + for j in range(ao_num_tot): + print(f'{i:5d} {j:5d} {i:5d} {j:5d} {Wfull[i,j,i,j].real:25.15E} {Wfull[i,j,i,j].imag:25.15E}') + print(f'{i:5d} {j:5d} {j:5d} {i:5d} {Wfull[i,j,j,i].real:25.15E} {Wfull[i,j,j,i].imag:25.15E}') for imo, iocc in enumerate(mo_occ): if iocc: E1 += 2*H1[imo,imo] diff --git a/src/ao_two_e_ints/cd_ao_ints.irp.f b/src/ao_two_e_ints/cd_ao_ints.irp.f index e1888bf9..73fee7c7 100644 --- a/src/ao_two_e_ints/cd_ao_ints.irp.f +++ b/src/ao_two_e_ints/cd_ao_ints.irp.f @@ -91,14 +91,14 @@ subroutine ao_map_fill_from_chol call idx2_tri_int(kj,kl,kjkl2) !TODO: verify the kj, kl as 4th index in expressions below if (kpt_sparse_map(kQ) > 0) then - ints_jl = chol_ao_integrals_complex(:,:,:,kl,kpt_sparse_map(kQ)) - !ints_jl = dconjg(chol_ao_integrals_complex(:,:,:,kl,kpt_sparse_map(kQ))) + !ints_jl = chol_ao_integrals_complex(:,:,:,kl,kpt_sparse_map(kQ)) + ints_jl = dconjg(chol_ao_integrals_complex(:,:,:,kl,kpt_sparse_map(kQ))) else do i_ao=1,ao_num_per_kpt do j_ao=1,ao_num_per_kpt do i_cd=1,chol_num_max - ints_jl(i_ao,j_ao,i_cd) = dconjg(chol_ao_integrals_complex(j_ao,i_ao,i_cd,kj,-kpt_sparse_map(kQ))) - !ints_jl(i_ao,j_ao,i_cd) = chol_ao_integrals_complex(j_ao,i_ao,i_cd,kj,-kpt_sparse_map(kQ)) + !ints_jl(i_ao,j_ao,i_cd) = dconjg(chol_ao_integrals_complex(j_ao,i_ao,i_cd,kj,-kpt_sparse_map(kQ))) + ints_jl(i_ao,j_ao,i_cd) = chol_ao_integrals_complex(j_ao,i_ao,i_cd,kj,-kpt_sparse_map(kQ)) enddo enddo enddo diff --git a/src/mo_two_e_ints/cd_mo_ints.irp.f b/src/mo_two_e_ints/cd_mo_ints.irp.f index 1afc2729..a832d510 100644 --- a/src/mo_two_e_ints/cd_mo_ints.irp.f +++ b/src/mo_two_e_ints/cd_mo_ints.irp.f @@ -160,7 +160,8 @@ subroutine mo_map_fill_from_chol_dot call idx2_tri_int(i,k,ik2) if (ik2 > jl2) exit !integral = zdotc(df_num,ints_jl(1,ij,il),1,ints_ik(1,ii,ik),1) - integral = zdotu(chol_num(kQ),ints_jl(1,ij,il),1,ints_ik(1,ii,ik),1) + !integral = zdotu(chol_num(kQ),ints_jl(1,ij,il),1,ints_ik(1,ii,ik),1) + integral = zdotu(chol_num(kQ),ints_jl(1,il,ij),1,ints_ik(1,ii,ik),1) ! print*,i,k,j,l,real(integral),imag(integral) if (cdabs(integral) < mo_integrals_threshold) then cycle diff --git a/src/utils_complex/dump_ao_2e_cplx.irp.f b/src/utils_complex/dump_ao_2e_cplx.irp.f index 2db5f614..63d81914 100644 --- a/src/utils_complex/dump_ao_2e_cplx.irp.f +++ b/src/utils_complex/dump_ao_2e_cplx.irp.f @@ -15,7 +15,9 @@ subroutine run do k=1,ao_num do l=1,ao_num tmp_cmplx = get_ao_two_e_integral_complex(i,j,k,l,ao_integrals_map,ao_integrals_map_2) - print'(4(I4),2(E23.15))',i,j,k,l,tmp_cmplx + if (cdabs(tmp_cmplx) .gt. 1E-10) then + print'(4(I4),2(E23.15))',i,j,k,l,tmp_cmplx + endif enddo enddo enddo From 4b4235d1619fe2260a6730db30053140c1148e9c Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Wed, 28 Sep 2022 10:54:01 -0500 Subject: [PATCH 25/28] separate jk terms for debugging --- src/hartree_fock/hf_energy.irp.f | 19 +- src/hartree_fock/print_e_scf.irp.f | 3 + src/scf_utils/fock_matrix_cplx.irp.f | 293 +++++++++++++++++++++++++-- 3 files changed, 296 insertions(+), 19 deletions(-) diff --git a/src/hartree_fock/hf_energy.irp.f b/src/hartree_fock/hf_energy.irp.f index 5a68164f..ea5e4b1c 100644 --- a/src/hartree_fock/hf_energy.irp.f +++ b/src/hartree_fock/hf_energy.irp.f @@ -13,19 +13,22 @@ END_PROVIDER BEGIN_PROVIDER [ double precision, hf_energy] &BEGIN_PROVIDER [ double precision, hf_two_electron_energy] +&BEGIN_PROVIDER [ double precision, hf_two_electron_energy_jk, (2)] &BEGIN_PROVIDER [ double precision, hf_one_electron_energy] implicit none BEGIN_DOC ! Hartree-Fock energy containing the nuclear repulsion, and its one- and two-body components. END_DOC - integer :: i,j,k + integer :: i,j,k,jk hf_energy = nuclear_repulsion hf_two_electron_energy = 0.d0 + hf_two_electron_energy_jk = 0.d0 hf_one_electron_energy = 0.d0 if (is_complex) then - complex*16 :: hf_1e_tmp, hf_2e_tmp + complex*16 :: hf_1e_tmp, hf_2e_tmp, hf_2e_tmp_jk(2) hf_1e_tmp = (0.d0,0.d0) hf_2e_tmp = (0.d0,0.d0) + hf_2e_tmp_jk = (0.d0,0.d0) do k=1,kpt_num do j=1,ao_num_per_kpt do i=1,ao_num_per_kpt @@ -33,9 +36,21 @@ END_PROVIDER +ao_two_e_integral_beta_kpts(i,j,k) * scf_density_matrix_ao_beta_kpts(j,i,k) ) hf_1e_tmp += ao_one_e_integrals_kpts(i,j,k) * (scf_density_matrix_ao_alpha_kpts(j,i,k) & + scf_density_matrix_ao_beta_kpts (j,i,k) ) + do jk=1,2 + hf_2e_tmp_jk(jk) += 0.5d0 * ( ao_two_e_integral_alpha_kpts_jk(i,j,k,jk) * scf_density_matrix_ao_alpha_kpts(j,i,k) & + +ao_two_e_integral_beta_kpts_jk(i,j,k,jk) * scf_density_matrix_ao_beta_kpts(j,i,k) ) + enddo enddo enddo enddo + do jk=1,2 + if (dabs(dimag(hf_2e_tmp_jk(jk))).gt.1.d-10) then + print*,'HF_2e energy (jk) should be real:',jk,irp_here + stop -1 + else + hf_two_electron_energy_jk(jk) = dble(hf_2e_tmp_jk(jk)) + endif + enddo if (dabs(dimag(hf_2e_tmp)).gt.1.d-10) then print*,'HF_2e energy should be real:',irp_here stop -1 diff --git a/src/hartree_fock/print_e_scf.irp.f b/src/hartree_fock/print_e_scf.irp.f index 989c0b9c..8efeba0d 100644 --- a/src/hartree_fock/print_e_scf.irp.f +++ b/src/hartree_fock/print_e_scf.irp.f @@ -15,6 +15,9 @@ subroutine run print*,hf_one_electron_energy print*,hf_two_electron_energy print*,hf_energy + print*,'hf 2e J,K energy' + print*,hf_two_electron_energy_jk(1) + print*,hf_two_electron_energy_jk(2) end diff --git a/src/scf_utils/fock_matrix_cplx.irp.f b/src/scf_utils/fock_matrix_cplx.irp.f index e2ada6fc..2f7d4a1a 100644 --- a/src/scf_utils/fock_matrix_cplx.irp.f +++ b/src/scf_utils/fock_matrix_cplx.irp.f @@ -542,27 +542,26 @@ BEGIN_PROVIDER [ complex*16, Fock_matrix_ao_kpts, (ao_num_per_kpt, ao_num_per_kp endif END_PROVIDER - - BEGIN_PROVIDER [ complex*16, ao_two_e_integral_alpha_kpts, (ao_num_per_kpt, ao_num_per_kpt, kpt_num) ] -&BEGIN_PROVIDER [ complex*16, ao_two_e_integral_beta_kpts , (ao_num_per_kpt, ao_num_per_kpt, kpt_num) ] + BEGIN_PROVIDER [ complex*16, ao_two_e_integral_alpha_kpts_jk, (ao_num_per_kpt, ao_num_per_kpt, kpt_num, 2) ] +&BEGIN_PROVIDER [ complex*16, ao_two_e_integral_beta_kpts_jk , (ao_num_per_kpt, ao_num_per_kpt, kpt_num, 2) ] use map_module implicit none BEGIN_DOC - ! Alpha and Beta Fock matrices in AO basis set + ! Alpha and Beta Fock matrices in AO basis set separated into j/k END_DOC !TODO: finish implementing this: see complex qp1 (different mapping) - + integer :: i,j,k,l,k1,r,s integer :: i0,j0,k0,l0 integer*8 :: p,q complex*16 :: integral, c0 - complex*16, allocatable :: ao_two_e_integral_alpha_tmp(:,:,:) - complex*16, allocatable :: ao_two_e_integral_beta_tmp(:,:,:) - - ao_two_e_integral_alpha_kpts = (0.d0,0.d0) - ao_two_e_integral_beta_kpts = (0.d0,0.d0) + complex*16, allocatable :: ao_two_e_integral_alpha_tmp(:,:,:,:) + complex*16, allocatable :: ao_two_e_integral_beta_tmp(:,:,:,:) + + ao_two_e_integral_alpha_kpts_jk = (0.d0,0.d0) + ao_two_e_integral_beta_kpts_jk = (0.d0,0.d0) PROVIDE ao_two_e_integrals_in_map scf_density_matrix_ao_alpha_kpts scf_density_matrix_ao_beta_kpts - + integer(omp_lock_kind) :: lck(ao_num) integer(map_size_kind) :: i8 integer :: ii(4), jj(4), kk(4), ll(4), k2 @@ -572,7 +571,267 @@ END_PROVIDER complex*16, parameter :: i_sign(4) = (/(0.d0,1.d0),(0.d0,1.d0),(0.d0,-1.d0),(0.d0,-1.d0)/) integer(key_kind) :: key1 integer :: kpt_i,kpt_j,kpt_k,kpt_l,idx_i,idx_j,idx_k,idx_l - + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,n_elements_max, & + !$OMP n_elements,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp, & + !$OMP kpt_i,kpt_j,kpt_k,kpt_l,idx_i,idx_j,idx_k,idx_l, & + !$OMP c0,key1)& + !$OMP SHARED(ao_num_per_kpt,SCF_density_matrix_ao_alpha_kpts, kpt_num, irp_here, & + !$OMP SCF_density_matrix_ao_beta_kpts, & + !$OMP ao_integrals_map, ao_two_e_integral_alpha_kpts_jk, ao_two_e_integral_beta_kpts_jk) + + call get_cache_map_n_elements_max(ao_integrals_map,n_elements_max) + allocate(keys(n_elements_max), values(n_elements_max)) + allocate(ao_two_e_integral_alpha_tmp(ao_num_per_kpt,ao_num_per_kpt,kpt_num,2), & + ao_two_e_integral_beta_tmp(ao_num_per_kpt,ao_num_per_kpt,kpt_num,2)) + ao_two_e_integral_alpha_tmp = (0.d0,0.d0) + ao_two_e_integral_beta_tmp = (0.d0,0.d0) + + !$OMP DO SCHEDULE(static,1) + do i8=0_8,ao_integrals_map%map_size + n_elements = n_elements_max + call get_cache_map(ao_integrals_map,i8,keys,values,n_elements) + do k1=1,n_elements + ! get original key + ! reverse of 2*key (imag part) and 2*key-1 (real part) + key1 = shiftr(keys(k1)+1,1) + + call two_e_integrals_index_reverse_complex_1(ii,jj,kk,ll,key1) + ! i<=k, j<=l, ik<=jl + ! ijkl, jilk, klij*, lkji* + + if (shiftl(key1,1)==keys(k1)) then !imaginary part (even) + do k2=1,4 + if (ii(k2)==0) then + cycle + endif + i = ii(k2) + j = jj(k2) + k = kk(k2) + l = ll(k2) + call get_kpt_idx_ao(i,kpt_i,idx_i) + call get_kpt_idx_ao(j,kpt_j,idx_j) + call get_kpt_idx_ao(k,kpt_k,idx_k) + call get_kpt_idx_ao(l,kpt_l,idx_l) + integral = i_sign(k2)*values(k1) !for klij and lkji, take complex conjugate + + !G_a(i,k) += D_{ab}(l,j)*() + !G_b(i,k) += D_{ab}(l,j)*() + !G_a(i,l) -= D_a (k,j)*() + !G_b(i,l) -= D_b (k,j)*() + + if (kpt_l.eq.kpt_j) then + c0 = (scf_density_matrix_ao_alpha_kpts(idx_l,idx_j,kpt_j)+scf_density_matrix_ao_beta_kpts(idx_l,idx_j,kpt_j))*integral + if(kpt_i.ne.kpt_k) then + print*,'problem in ',irp_here,' ikjl: ',kpt_i,kpt_k,kpt_j,kpt_l + stop 1 + endif + ao_two_e_integral_alpha_tmp(idx_i,idx_k,kpt_i,1) += c0 + ao_two_e_integral_beta_tmp (idx_i,idx_k,kpt_i,1) += c0 + endif + + if (kpt_l.eq.kpt_i) then + if(kpt_j.ne.kpt_k) then + print*,'problem in ',irp_here,' ikjl: ',kpt_i,kpt_k,kpt_j,kpt_l + stop 1 + endif + ao_two_e_integral_alpha_tmp(idx_i,idx_l,kpt_i,2) -= SCF_density_matrix_ao_alpha_kpts(idx_k,idx_j,kpt_j) * integral + ao_two_e_integral_beta_tmp (idx_i,idx_l,kpt_i,2) -= scf_density_matrix_ao_beta_kpts (idx_k,idx_j,kpt_j) * integral + endif + enddo + else ! real part + do k2=1,4 + if (ii(k2)==0) then + cycle + endif + i = ii(k2) + j = jj(k2) + k = kk(k2) + l = ll(k2) + call get_kpt_idx_ao(i,kpt_i,idx_i) + call get_kpt_idx_ao(j,kpt_j,idx_j) + call get_kpt_idx_ao(k,kpt_k,idx_k) + call get_kpt_idx_ao(l,kpt_l,idx_l) + integral = values(k1) + + if (kpt_l.eq.kpt_j) then + c0 = (scf_density_matrix_ao_alpha_kpts(idx_l,idx_j,kpt_j)+scf_density_matrix_ao_beta_kpts(idx_l,idx_j,kpt_j))*integral + if(kpt_i.ne.kpt_k) then + print*,'problem in ',irp_here,' ikjl: ',kpt_i,kpt_k,kpt_j,kpt_l + stop 1 + endif + ao_two_e_integral_alpha_tmp(idx_i,idx_k,kpt_i,1) += c0 + ao_two_e_integral_beta_tmp (idx_i,idx_k,kpt_i,1) += c0 + endif + + if (kpt_l.eq.kpt_i) then + if(kpt_j.ne.kpt_k) then + print*,'problem in ',irp_here,' ikjl: ',kpt_i,kpt_k,kpt_j,kpt_l + stop 1 + endif + ao_two_e_integral_alpha_tmp(idx_i,idx_l,kpt_i,2) -= SCF_density_matrix_ao_alpha_kpts(idx_k,idx_j,kpt_j) * integral + ao_two_e_integral_beta_tmp (idx_i,idx_l,kpt_i,2) -= scf_density_matrix_ao_beta_kpts (idx_k,idx_j,kpt_j) * integral + endif + enddo + endif + enddo + enddo + !$OMP END DO NOWAIT + !$OMP CRITICAL + ao_two_e_integral_alpha_kpts_jk += ao_two_e_integral_alpha_tmp + ao_two_e_integral_beta_kpts_jk += ao_two_e_integral_beta_tmp + !$OMP END CRITICAL + deallocate(keys,values,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp) + !$OMP END PARALLEL + + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,n_elements_max, & + !$OMP n_elements,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp, & + !$OMP kpt_i,kpt_j,kpt_k,kpt_l,idx_i,idx_j,idx_k,idx_l, & + !$OMP c0,key1)& + !$OMP SHARED(ao_num_per_kpt,SCF_density_matrix_ao_alpha_kpts,kpt_num, irp_here, & + !$OMP SCF_density_matrix_ao_beta_kpts, & + !$OMP ao_integrals_map_2, ao_two_e_integral_alpha_kpts_jk, ao_two_e_integral_beta_kpts_jk) + + call get_cache_map_n_elements_max(ao_integrals_map_2,n_elements_max) + allocate(keys(n_elements_max), values(n_elements_max)) + allocate(ao_two_e_integral_alpha_tmp(ao_num_per_kpt,ao_num_per_kpt,kpt_num,2), & + ao_two_e_integral_beta_tmp(ao_num_per_kpt,ao_num_per_kpt,kpt_num,2)) + ao_two_e_integral_alpha_tmp = (0.d0,0.d0) + ao_two_e_integral_beta_tmp = (0.d0,0.d0) + + !$OMP DO SCHEDULE(static,1) + do i8=0_8,ao_integrals_map_2%map_size + n_elements = n_elements_max + call get_cache_map(ao_integrals_map_2,i8,keys,values,n_elements) + do k1=1,n_elements + ! get original key + ! reverse of 2*key (imag part) and 2*key-1 (real part) + key1 = shiftr(keys(k1)+1,1) + + call two_e_integrals_index_reverse_complex_2(ii,jj,kk,ll,key1) + ! i>=k, j<=l, ik<=jl + ! ijkl, jilk, klij*, lkji* + if (shiftl(key1,1)==keys(k1)) then !imaginary part + do k2=1,4 + if (ii(k2)==0) then + cycle + endif + i = ii(k2) + j = jj(k2) + k = kk(k2) + l = ll(k2) + call get_kpt_idx_ao(i,kpt_i,idx_i) + call get_kpt_idx_ao(j,kpt_j,idx_j) + call get_kpt_idx_ao(k,kpt_k,idx_k) + call get_kpt_idx_ao(l,kpt_l,idx_l) + integral = i_sign(k2)*values(k1) ! for klij and lkji, take conjugate + + !G_a(i,k) += D_{ab}(l,j)*() + !G_b(i,k) += D_{ab}(l,j)*() + !G_a(i,l) -= D_a (k,j)*() + !G_b(i,l) -= D_b (k,j)*() + + if (kpt_l.eq.kpt_j) then + c0 = (scf_density_matrix_ao_alpha_kpts(idx_l,idx_j,kpt_j)+scf_density_matrix_ao_beta_kpts(idx_l,idx_j,kpt_j))*integral + if(kpt_i.ne.kpt_k) then + print*,'problem in ',irp_here,' ikjl: ',kpt_i,kpt_k,kpt_j,kpt_l + stop 1 + endif + ao_two_e_integral_alpha_tmp(idx_i,idx_k,kpt_i,1) += c0 + ao_two_e_integral_beta_tmp (idx_i,idx_k,kpt_i,1) += c0 + endif + + if (kpt_l.eq.kpt_i) then + if(kpt_j.ne.kpt_k) then + print*,'problem in ',irp_here,' ikjl: ',kpt_i,kpt_k,kpt_j,kpt_l + stop 1 + endif + ao_two_e_integral_alpha_tmp(idx_i,idx_l,kpt_i,2) -= SCF_density_matrix_ao_alpha_kpts(idx_k,idx_j,kpt_j) * integral + ao_two_e_integral_beta_tmp (idx_i,idx_l,kpt_i,2) -= scf_density_matrix_ao_beta_kpts (idx_k,idx_j,kpt_j) * integral + endif + enddo + else ! real part + do k2=1,4 + if (ii(k2)==0) then + cycle + endif + i = ii(k2) + j = jj(k2) + k = kk(k2) + l = ll(k2) + call get_kpt_idx_ao(i,kpt_i,idx_i) + call get_kpt_idx_ao(j,kpt_j,idx_j) + call get_kpt_idx_ao(k,kpt_k,idx_k) + call get_kpt_idx_ao(l,kpt_l,idx_l) + integral = values(k1) + + if (kpt_l.eq.kpt_j) then + c0 = (scf_density_matrix_ao_alpha_kpts(idx_l,idx_j,kpt_j)+scf_density_matrix_ao_beta_kpts(idx_l,idx_j,kpt_j))*integral + if(kpt_i.ne.kpt_k) then + print*,'problem in ',irp_here,' ikjl: ',kpt_i,kpt_k,kpt_j,kpt_l + stop 1 + endif + ao_two_e_integral_alpha_tmp(idx_i,idx_k,kpt_i,1) += c0 + ao_two_e_integral_beta_tmp (idx_i,idx_k,kpt_i,1) += c0 + endif + + if (kpt_l.eq.kpt_i) then + if(kpt_j.ne.kpt_k) then + print*,'problem in ',irp_here,' ikjl: ',kpt_i,kpt_k,kpt_j,kpt_l + stop 1 + endif + ao_two_e_integral_alpha_tmp(idx_i,idx_l,kpt_i,2) -= SCF_density_matrix_ao_alpha_kpts(idx_k,idx_j,kpt_j) * integral + ao_two_e_integral_beta_tmp (idx_i,idx_l,kpt_i,2) -= scf_density_matrix_ao_beta_kpts (idx_k,idx_j,kpt_j) * integral + endif + enddo + endif + enddo + enddo + !$OMP END DO NOWAIT + !$OMP CRITICAL + ao_two_e_integral_alpha_kpts_jk += ao_two_e_integral_alpha_tmp + ao_two_e_integral_beta_kpts_jk += ao_two_e_integral_beta_tmp + !$OMP END CRITICAL + deallocate(keys,values,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp) + !$OMP END PARALLEL + + +END_PROVIDER + + + BEGIN_PROVIDER [ complex*16, ao_two_e_integral_alpha_kpts, (ao_num_per_kpt, ao_num_per_kpt, kpt_num) ] +&BEGIN_PROVIDER [ complex*16, ao_two_e_integral_beta_kpts , (ao_num_per_kpt, ao_num_per_kpt, kpt_num) ] + use map_module + implicit none + BEGIN_DOC + ! Alpha and Beta Fock matrices in AO basis set + END_DOC + !TODO: finish implementing this: see complex qp1 (different mapping) + + integer :: i,j,k,l,k1,r,s + integer :: i0,j0,k0,l0 + integer*8 :: p,q + complex*16 :: integral, c0 + complex*16, allocatable :: ao_two_e_integral_alpha_tmp(:,:,:) + complex*16, allocatable :: ao_two_e_integral_beta_tmp(:,:,:) + + ao_two_e_integral_alpha_kpts = (0.d0,0.d0) + ao_two_e_integral_beta_kpts = (0.d0,0.d0) + PROVIDE ao_two_e_integrals_in_map scf_density_matrix_ao_alpha_kpts scf_density_matrix_ao_beta_kpts + + integer(omp_lock_kind) :: lck(ao_num) + integer(map_size_kind) :: i8 + integer :: ii(4), jj(4), kk(4), ll(4), k2 + integer(cache_map_size_kind) :: n_elements_max, n_elements + integer(key_kind), allocatable :: keys(:) + double precision, allocatable :: values(:) + complex*16, parameter :: i_sign(4) = (/(0.d0,1.d0),(0.d0,1.d0),(0.d0,-1.d0),(0.d0,-1.d0)/) + integer(key_kind) :: key1 + integer :: kpt_i,kpt_j,kpt_k,kpt_l,idx_i,idx_j,idx_k,idx_l + !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,n_elements_max, & !$OMP n_elements,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp, & @@ -581,14 +840,14 @@ END_PROVIDER !$OMP SHARED(ao_num_per_kpt,SCF_density_matrix_ao_alpha_kpts, kpt_num, irp_here, & !$OMP SCF_density_matrix_ao_beta_kpts, & !$OMP ao_integrals_map, ao_two_e_integral_alpha_kpts, ao_two_e_integral_beta_kpts) - + call get_cache_map_n_elements_max(ao_integrals_map,n_elements_max) allocate(keys(n_elements_max), values(n_elements_max)) allocate(ao_two_e_integral_alpha_tmp(ao_num_per_kpt,ao_num_per_kpt,kpt_num), & ao_two_e_integral_beta_tmp(ao_num_per_kpt,ao_num_per_kpt,kpt_num)) ao_two_e_integral_alpha_tmp = (0.d0,0.d0) ao_two_e_integral_beta_tmp = (0.d0,0.d0) - + !$OMP DO SCHEDULE(static,1) do i8=0_8,ao_integrals_map%map_size n_elements = n_elements_max @@ -686,7 +945,7 @@ END_PROVIDER deallocate(keys,values,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp) !$OMP END PARALLEL - + !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,n_elements_max, & !$OMP n_elements,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp, & @@ -695,14 +954,14 @@ END_PROVIDER !$OMP SHARED(ao_num_per_kpt,SCF_density_matrix_ao_alpha_kpts,kpt_num, irp_here, & !$OMP SCF_density_matrix_ao_beta_kpts, & !$OMP ao_integrals_map_2, ao_two_e_integral_alpha_kpts, ao_two_e_integral_beta_kpts) - + call get_cache_map_n_elements_max(ao_integrals_map_2,n_elements_max) allocate(keys(n_elements_max), values(n_elements_max)) allocate(ao_two_e_integral_alpha_tmp(ao_num_per_kpt,ao_num_per_kpt,kpt_num), & ao_two_e_integral_beta_tmp(ao_num_per_kpt,ao_num_per_kpt,kpt_num)) ao_two_e_integral_alpha_tmp = (0.d0,0.d0) ao_two_e_integral_beta_tmp = (0.d0,0.d0) - + !$OMP DO SCHEDULE(static,1) do i8=0_8,ao_integrals_map_2%map_size n_elements = n_elements_max From 7486367048d0031e148e6a77a8d9e63c8c37086b Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Wed, 28 Sep 2022 11:21:01 -0500 Subject: [PATCH 26/28] change to mulliken ordering when dumping ao ints --- src/utils_complex/dump_ao_2e_cplx.irp.f | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/utils_complex/dump_ao_2e_cplx.irp.f b/src/utils_complex/dump_ao_2e_cplx.irp.f index 63d81914..ea36b876 100644 --- a/src/utils_complex/dump_ao_2e_cplx.irp.f +++ b/src/utils_complex/dump_ao_2e_cplx.irp.f @@ -16,20 +16,20 @@ subroutine run do l=1,ao_num tmp_cmplx = get_ao_two_e_integral_complex(i,j,k,l,ao_integrals_map,ao_integrals_map_2) if (cdabs(tmp_cmplx) .gt. 1E-10) then - print'(4(I4),2(E23.15))',i,j,k,l,tmp_cmplx + print'(4(I4),2(E23.15))',i,k,j,l,tmp_cmplx endif enddo enddo enddo enddo - print*,'map1' - do i=0,ao_integrals_map%map_size - print*,i,ao_integrals_map%map(i)%value(:) - print*,i,ao_integrals_map%map(i)%key(:) - enddo - print*,'map2' - do i=0,ao_integrals_map_2%map_size - print*,i,ao_integrals_map_2%map(i)%value(:) - print*,i,ao_integrals_map_2%map(i)%key(:) - enddo + !print*,'map1' + !do i=0,ao_integrals_map%map_size + ! print*,i,ao_integrals_map%map(i)%value(:) + ! print*,i,ao_integrals_map%map(i)%key(:) + !enddo + !print*,'map2' + !do i=0,ao_integrals_map_2%map_size + ! print*,i,ao_integrals_map_2%map(i)%value(:) + ! print*,i,ao_integrals_map_2%map(i)%key(:) + !enddo end From 914ef54417d4ecc85f52dbfdb7873c1d5710c178 Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Wed, 5 Oct 2022 15:54:29 -0500 Subject: [PATCH 27/28] idx perm fix --- bin/qp_convert_h5_to_ezfio | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/bin/qp_convert_h5_to_ezfio b/bin/qp_convert_h5_to_ezfio index 4bee742f..ad5c4dbe 100755 --- a/bin/qp_convert_h5_to_ezfio +++ b/bin/qp_convert_h5_to_ezfio @@ -613,14 +613,23 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True,is_ao=True): # B = LQ0[ks].conj() # W = np.einsum('prn,sqn->pqrs',A,B) # 5 + #if Qneg: + # A = LQ0[kr].transpose((1,0,2)).conj() + # B = LQ0[kq].transpose((1,0,2)) + # W = np.einsum('prn,qsn->pqrs',A,B) + #else: + # A = LQ0[kp] + # B = LQ0[ks].conj() + # W = np.einsum('prn,qsn->pqrs',A,B) + # 6 if Qneg: A = LQ0[kr].transpose((1,0,2)).conj() B = LQ0[kq].transpose((1,0,2)) - W = np.einsum('prn,qsn->pqrs',A,B) + W = np.einsum('prn,sqn->pqrs',A,B) else: A = LQ0[kp] B = LQ0[ks].conj() - W = np.einsum('prn,qsn->pqrs',A,B) + W = np.einsum('prn,sqn->pqrs',A,B) p0 = kp*ao_num_per_kpt r0 = kr*ao_num_per_kpt q0 = kq*ao_num_per_kpt @@ -630,7 +639,8 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True,is_ao=True): for ir in range(ao_num_per_kpt): for i_s in range(ao_num_per_kpt): v = W[ip,iq,ir,i_s] - print(f'{p0+ip:5d} {q0+iq:5d} {r0+ir:5d} {s0+i_s:5d} {v.real:25.15E} {v.imag:25.15E}') + #print(f'{p0+ip:5d} {q0+iq:5d} {r0+ir:5d} {s0+i_s:5d} {v.real:25.15E} {v.imag:25.15E}') + print(f'{p0+ip:5d} {r0+ir:5d} {q0+iq:5d} {s0+i_s:5d} {v.real:25.15E} {v.imag:25.15E}') Wfull[p0:p0+ao_num_per_kpt,q0:q0+ao_num_per_kpt,r0:r0+ao_num_per_kpt,s0:s0+ao_num_per_kpt] = W.copy() H1 = np.zeros((ao_num_tot, ao_num_tot), dtype=np.complex128) for Qi in range(kpt_num): From 199d5398f46e82beec33f66bbe15af00442555f7 Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Thu, 6 Oct 2022 12:23:59 -0500 Subject: [PATCH 28/28] cleanup --- bin/qp_convert_h5_to_ezfio | 19 ++++++++----------- 1 file changed, 8 insertions(+), 11 deletions(-) diff --git a/bin/qp_convert_h5_to_ezfio b/bin/qp_convert_h5_to_ezfio index ad5c4dbe..f3e78d53 100755 --- a/bin/qp_convert_h5_to_ezfio +++ b/bin/qp_convert_h5_to_ezfio @@ -288,11 +288,7 @@ def convert_mol(filename,qph5path): def convert_kpts_cd(filename,qph5path,qmcpack=True,is_ao=True): import json from scipy.linalg import block_diag - #dump_fci, dump_cd = (False, False) - dump_fci, dump_cd = (True, False) - #dump_fci, dump_cd = (False, True) - #dump_fci2 = True - dump_fci2 = False + dump_fci, dump_cd, dump_fci2 = (False, False, False) ezfio.set_file(filename) ezfio.set_nuclei_is_complex(True) @@ -495,11 +491,11 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True,is_ao=True): # should this be in ao_basis? ao_two_e_ints? with h5py.File(qph5path,'r') as qph5: nchol_per_kpt_all = qph5['Hamiltonian']['NCholPerKP'][:] - print(nchol_per_kpt_all) + print(f'nchol_per_kpt_full = {nchol_per_kpt_all}') #nchol_per_kpt = nchol_per_kpt_all[nchol_per_kpt_all != 0] nchol_per_kpt = nchol_per_kpt_all[np.array(unique_k_idx,dtype=int)-1] - print(nchol_per_kpt) - print(unique_k_idx) + print(f'nchol_per_kpt = {nchol_per_kpt}') + print(f'unique_k_idx = {unique_k_idx}') #for i in range(kpt_num): # if i+1 in unique_k_idx: # print('* ',i,nchol_per_kpt_all[i]) @@ -515,11 +511,11 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True,is_ao=True): #ao_chol_two_e_ints = np.zeros((2, ao_num_per_kpt, ao_num_per_kpt, nchol_per_kpt_max, kpt_num, len(nchol_per_kpt))) L_list = [] L_all = np.zeros((unique_kpt_num, kpt_num, ao_num_per_kpt, ao_num_per_kpt, nchol_per_kpt_max,2),dtype=np.float64) - print(kpt_sparse_map) - print(np.array(unique_k_idx)-1) + print(f'kpt_sparse_map = {kpt_sparse_map}') + print(f'unique_k_idx-1 = {np.array(unique_k_idx)-1}') for i in range(unique_kpt_num): ki = unique_k_idx[i]-1 - print(i, ki) + #print(i, ki) L_i = qph5[f'Hamiltonian/KPFactorized/L{ki}'][()].reshape((kpt_num, ao_num_per_kpt, ao_num_per_kpt, nchol_per_kpt[kpt_sparse_map[ki]-1], 2)) #L.reshape(kpt_num, ao_num_per_kpt, ao_num_per_kpt, nchol_per_kpt[i], 2) #L = np.einsum("ijklm->ilkjm", A, B) @@ -839,6 +835,7 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True,is_ao=True): ezfio.set_ao_two_e_ints_io_df_ao_integrals('Read') """ else: + raise NotImplementedError """ ezfio.set_io_chol_mo_integrals('Read') df_num = qph5['ao_two_e_ints'].attrs['df_num']