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