From 1daadf2dcfbd24d6ff7c0e9688aac1741acbb196 Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Tue, 2 Jun 2020 11:55:33 -0500 Subject: [PATCH 1/8] added gf module; initial commit --- src/green/EZFIO.cfg | 101 ++++ src/green/NEED | 1 + src/green/README.rst | 6 + src/green/green.main.irp.f | 51 ++ src/green/hu0_hp.irp.f | 847 ++++++++++++++++++++++++++ src/green/hu0_lanczos.irp.f | 405 +++++++++++++ src/green/lanczos.irp.f | 882 ++++++++++++++++++++++++++++ src/green/plot-spec-dens.py | 90 +++ src/green/print_dets_test.irp.f | 15 + src/green/print_e_mo_debug.irp.f | 15 + src/green/print_h_debug.irp.f | 178 ++++++ src/green/print_h_omp_debug.irp.f | 41 ++ src/green/print_spectral_dens.irp.f | 43 ++ src/green/utils_hp.irp.f | 614 +++++++++++++++++++ 14 files changed, 3289 insertions(+) create mode 100644 src/green/EZFIO.cfg create mode 100644 src/green/NEED create mode 100644 src/green/README.rst create mode 100644 src/green/green.main.irp.f create mode 100644 src/green/hu0_hp.irp.f create mode 100644 src/green/hu0_lanczos.irp.f create mode 100644 src/green/lanczos.irp.f create mode 100755 src/green/plot-spec-dens.py create mode 100644 src/green/print_dets_test.irp.f create mode 100644 src/green/print_e_mo_debug.irp.f create mode 100644 src/green/print_h_debug.irp.f create mode 100644 src/green/print_h_omp_debug.irp.f create mode 100644 src/green/print_spectral_dens.irp.f create mode 100644 src/green/utils_hp.irp.f diff --git a/src/green/EZFIO.cfg b/src/green/EZFIO.cfg new file mode 100644 index 00000000..859a3eb9 --- /dev/null +++ b/src/green/EZFIO.cfg @@ -0,0 +1,101 @@ +[n_lanczos_complete] +type: integer +doc: number of lanczos iterations completed +interface: ezfio,provider,ocaml +default: 0 + +[n_lanczos_iter] +type: integer +doc: number of lanczos iterations +interface: ezfio,provider,ocaml +default: 10 + +[omega_min] +type: double precision +doc: lower limit of frequency for spectral density calculation +interface: ezfio,provider,ocaml +default: -2.e-1 + +[omega_max] +type: double precision +doc: upper limit of frequency for spectral density calculation +interface: ezfio,provider,ocaml +default: 1.2e1 + +[n_omega] +type: integer +doc: number of points for spectral density calculation +interface: ezfio,provider,ocaml +default: 1000 + +[gf_epsilon] +type: double precision +doc: infinitesimal imaginary frequency term in Green's function +interface: ezfio,provider,ocaml +default: 1.e-2 + +[n_green_vec] +type: integer +doc: number of holes/particles +interface: ezfio +default: 2 + +[green_idx] +interface: ezfio +doc: homo/lumo indices +type: integer +size: (green.n_green_vec) + +[green_spin] +interface: ezfio +doc: homo/lumo spin +type: integer +size: (green.n_green_vec) + +[green_sign] +interface: ezfio +doc: homo/lumo sign +type: double precision +size: (green.n_green_vec) + +[alpha_lanczos] +interface: ezfio +doc: lanczos alpha values +type: double precision +size: (green.n_lanczos_iter,green.n_green_vec) + +[beta_lanczos] +interface: ezfio +doc: lanczos beta values +type: double precision +size: (green.n_lanczos_iter,green.n_green_vec) + +[un_lanczos] +interface: ezfio +doc: saved lanczos u vector +type: complex*16 +size: (determinants.n_det,green.n_green_vec) + +[vn_lanczos] +interface: ezfio +doc: saved lanczos v vector +type: complex*16 +size: (determinants.n_det,green.n_green_vec) + +[lanczos_eigvals] +interface: ezfio +doc: eigvals of tridiagonal form of H +type: double precision +size: (green.n_lanczos_iter,green.n_green_vec) + +[lanczos_debug_print] +interface: ezfio,provider,ocaml +type: logical +doc: printing of lanczos vectors at every step +default: False + +[n_lanczos_debug] +interface: ezfio,provider,ocaml +type: integer +doc: number of elements to print +default: 10 diff --git a/src/green/NEED b/src/green/NEED new file mode 100644 index 00000000..4315d882 --- /dev/null +++ b/src/green/NEED @@ -0,0 +1 @@ +davidson fci diff --git a/src/green/README.rst b/src/green/README.rst new file mode 100644 index 00000000..6bdb2ca7 --- /dev/null +++ b/src/green/README.rst @@ -0,0 +1,6 @@ +===== +dummy +===== + +Module necessary to avoid the ``xxx is a root module but does not contain a main file`` message. + diff --git a/src/green/green.main.irp.f b/src/green/green.main.irp.f new file mode 100644 index 00000000..c9b3ef66 --- /dev/null +++ b/src/green/green.main.irp.f @@ -0,0 +1,51 @@ +program green + implicit none + BEGIN_DOC +! TODO + END_DOC + read_wf = .True. + touch read_wf + provide n_green_vec + print*,'ref_bitmask_energy = ',ref_bitmask_energy +! call psicoefprinttest + call print_lanczos_eigvals + call print_spec +end + +subroutine psicoefprinttest + implicit none + integer :: i + TOUCH psi_coef + print *, 'printing ndet', N_det +end +subroutine print_lanczos_eigvals + implicit none + integer :: i, iunit, j + integer :: getunitandopen + character(5) :: jstr + + do j=1,n_green_vec + write(jstr,'(I0.3)') j + iunit = getunitandopen('lanczos_eigval_alpha_beta.out.'//trim(jstr),'w') + print *, 'printing lanczos eigenvalues, alpha, beta to "lanczos_eigval_alpha_beta.out.'//trim(jstr)//'"' + do i=1,n_lanczos_iter + write(iunit,'(I6,3(E25.15))') i, lanczos_eigvals(i,j), alpha_lanczos(i,j), beta_lanczos(i,j) + enddo + close(iunit) + enddo +end +subroutine print_spec + implicit none + integer :: i, iunit, j + integer :: getunitandopen + character(5) :: jstr + do j=1,n_green_vec + write(jstr,'(I0.3)') j + iunit = getunitandopen('omega_A.out.'//trim(jstr),'w') + print *, 'printing frequency, spectral density to "omega_A.out.'//trim(jstr)//'"' + do i=1,n_omega + write(iunit,'(2(E25.15))') omega_list(i), spectral_lanczos(i,j) + enddo + close(iunit) + enddo +end diff --git a/src/green/hu0_hp.irp.f b/src/green/hu0_hp.irp.f new file mode 100644 index 00000000..c3d8be40 --- /dev/null +++ b/src/green/hu0_hp.irp.f @@ -0,0 +1,847 @@ +! modified from H_S2_u_0_nstates_openmp in Davidson/u0Hu0.irp.f + +subroutine h_u_0_hp_openmp(v_0,u_0,N_hp,sze,spin_hp,sign_hp,idx_hp) + use bitmasks + implicit none + BEGIN_DOC + ! Computes v_0 = H|u_0> + ! + ! Assumes that the determinants are in psi_det + ! + ! istart, iend, ishift, istep are used in ZMQ parallelization. + ! + ! N_hp is number of holes and particles to be applied + ! each element of spin_hp is either 1(alpha) or 2(beta) + ! each element of sign_hp is either 1(particle) or -1(hole) + ! idx_hp contains orbital indices for holes and particles + END_DOC + integer, intent(in) :: N_hp,sze + complex*16, intent(inout) :: v_0(sze,N_hp), u_0(sze,N_hp) + integer :: k + complex*16, allocatable :: u_t(:,:), v_t(:,:) + integer, intent(in) :: spin_hp(N_hp), idx_hp(N_hp) + double precision, intent(in) :: sign_hp(N_hp) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: u_t + allocate(u_t(N_hp,N_det),v_t(N_hp,N_det)) + do k=1,N_hp + call cdset_order(u_0(1,k),psi_bilinear_matrix_order,N_det) + enddo + v_t = (0.d0,0.d0) + call cdtranspose( & + u_0, & + size(u_0, 1), & + u_t, & + size(u_t, 1), & + N_det, N_hp) + + call h_u_0_hp_openmp_work(v_t,u_t,N_hp,sze,spin_hp,sign_hp,idx_hp,1,N_det,0,1) + deallocate(u_t) + + call cdtranspose( & + v_t, & + size(v_t, 1), & + v_0, & + size(v_0, 1), & + N_hp, N_det) + deallocate(v_t) + + do k=1,N_hp + call cdset_order(v_0(1,k),psi_bilinear_matrix_order_reverse,N_det) + call cdset_order(u_0(1,k),psi_bilinear_matrix_order_reverse,N_det) + enddo + +end + + +subroutine h_u_0_hp_openmp_work(v_t,u_t,N_hp,sze,spin_hp,sign_hp,idx_hp,istart,iend,ishift,istep) + use bitmasks + implicit none + BEGIN_DOC + ! Computes v_t = H|u_t> + ! + ! Default should be 1,N_det,0,1 + END_DOC + integer, intent(in) :: N_hp,sze,istart,iend,ishift,istep + complex*16, intent(in) :: u_t(N_hp,N_det) + complex*16, intent(out) :: v_t(N_hp,sze) + integer, intent(in) :: spin_hp(N_hp), idx_hp(N_hp) + double precision, intent(in) :: sign_hp(N_hp) + + + PROVIDE ref_bitmask_energy N_int + + select case (N_int) + case (1) + call H_u_0_hp_openmp_work_1(v_t,u_t,N_hp,sze,spin_hp,sign_hp,idx_hp,istart,iend,ishift,istep) + case (2) + call H_u_0_hp_openmp_work_2(v_t,u_t,N_hp,sze,spin_hp,sign_hp,idx_hp,istart,iend,ishift,istep) + case (3) + call H_u_0_hp_openmp_work_3(v_t,u_t,N_hp,sze,spin_hp,sign_hp,idx_hp,istart,iend,ishift,istep) + case (4) + call H_u_0_hp_openmp_work_4(v_t,u_t,N_hp,sze,spin_hp,sign_hp,idx_hp,istart,iend,ishift,istep) + case default + call H_u_0_hp_openmp_work_N_int(v_t,u_t,N_hp,sze,spin_hp,sign_hp,idx_hp,istart,iend,ishift,istep) + end select +end +BEGIN_TEMPLATE + +subroutine h_u_0_hp_openmp_work_$N_int(v_t,u_t,N_hp,sze,spin_hp,sign_hp,idx_hp,istart,iend,ishift,istep) + use bitmasks + implicit none + BEGIN_DOC + ! Computes v_t = H|u_t> and s_t = S^2 |u_t> + ! + ! Default should be 1,N_det,0,1 + END_DOC + integer, intent(in) :: N_hp,sze,istart,iend,ishift,istep + complex*16, intent(in) :: u_t(N_hp,N_det) + complex*16, intent(out) :: v_t(N_hp,sze) + integer, intent(in) :: spin_hp(N_hp), idx_hp(N_hp) + double precision, intent(in) :: sign_hp(N_hp) + + complex*16 :: hij + double precision :: hii + integer :: i,j,k,l + integer :: k_a, k_b, l_a, l_b, m_a, m_b + integer :: istate + integer :: krow, kcol, krow_b, kcol_b + integer :: lrow, lcol + integer :: mrow, mcol + integer(bit_kind) :: spindet($N_int) + integer(bit_kind) :: tmp_det($N_int,2) + integer(bit_kind) :: tmp_det2($N_int,2) + integer(bit_kind) :: tmp_det3($N_int,2) + integer(bit_kind), allocatable :: buffer(:,:) + integer :: n_doubles + integer, allocatable :: doubles(:) + integer, allocatable :: singles_a(:) + integer, allocatable :: singles_b(:) + integer, allocatable :: idx(:), idx0(:) + integer :: maxab, n_singles_a, n_singles_b, kcol_prev + integer*8 :: k8 + + logical, allocatable :: exc_is_banned_a1(:),exc_is_banned_b1(:),exc_is_banned_a2(:),exc_is_banned_b2(:) + logical, allocatable :: exc_is_banned_ab1(:),exc_is_banned_ab12(:),allowed_hp(:) + logical :: all_banned_a1,all_banned_b1,all_banned_a2,all_banned_b2 + logical :: all_banned_ab12,all_banned_ab1 + integer :: ii,na,nb + double precision, allocatable :: hii_hp(:) + complex*16, allocatable :: hij_hp(:) + + maxab = max(N_det_alpha_unique, N_det_beta_unique)+1 + allocate(idx0(maxab)) + + do i=1,maxab + idx0(i) = i + enddo + + ! Prepare the array of all alpha single excitations + ! ------------------------------------------------- + + PROVIDE N_int nthreads_davidson elec_num + !$OMP PARALLEL DEFAULT(NONE) NUM_THREADS(nthreads_davidson) & + !$OMP SHARED(psi_bilinear_matrix_rows, N_det, & + !$OMP psi_bilinear_matrix_columns, & + !$OMP psi_det_alpha_unique, psi_det_beta_unique, & + !$OMP n_det_alpha_unique, n_det_beta_unique, N_int, & + !$OMP psi_bilinear_matrix_transp_rows, & + !$OMP psi_bilinear_matrix_transp_columns, & + !$OMP psi_bilinear_matrix_transp_order, N_hp, & + !$OMP psi_bilinear_matrix_order_transp_reverse, & + !$OMP psi_bilinear_matrix_columns_loc, & + !$OMP psi_bilinear_matrix_transp_rows_loc, & + !$OMP istart, iend, istep, irp_here, v_t, & + !$OMP spin_hp,sign_hp,idx_hp, & + !$OMP elec_num_tab,nuclear_repulsion, & + !$OMP ishift, idx0, u_t, maxab) & + !$OMP PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i, & + !$OMP lcol, lrow, l_a, l_b, & + !$OMP buffer, doubles, n_doubles, & + !$OMP tmp_det2, hii, hij, idx, l, kcol_prev, & + !$OMP singles_a, n_singles_a, singles_b, & + !$OMP exc_is_banned_a1,exc_is_banned_b1,exc_is_banned_ab1, & + !$OMP exc_is_banned_a2,exc_is_banned_b2,exc_is_banned_ab12, & + !$OMP all_banned_a1,all_banned_b1,all_banned_ab1, & + !$OMP all_banned_a2,all_banned_b2,all_banned_ab12, & + !$OMP allowed_hp, & + !$OMP ii, hij_hp, j, hii_hp,na,nb, & + !$OMP n_singles_b, k8) + + ! Alpha/Beta double excitations + ! ============================= + + allocate( buffer($N_int,maxab), & + singles_a(maxab), & + singles_b(maxab), & + doubles(maxab), & + idx(maxab), & + exc_is_banned_a1(N_hp), & + exc_is_banned_b1(N_hp), & + exc_is_banned_a2(N_hp), & + exc_is_banned_b2(N_hp), & + exc_is_banned_ab1(N_hp), & + exc_is_banned_ab12(N_hp), & + allowed_hp(N_hp), & + hij_hp(N_hp), & + hii_hp(N_hp)) + + kcol_prev=-1 + all_banned_b1=.False. + ASSERT (iend <= N_det) + ASSERT (istart > 0) + ASSERT (istep > 0) + + !$OMP DO SCHEDULE(dynamic,64) + do k_a=istart+ishift,iend,istep + ! iterate over dets in psi + + krow = psi_bilinear_matrix_rows(k_a) + ASSERT (krow <= N_det_alpha_unique) + + kcol = psi_bilinear_matrix_columns(k_a) + ASSERT (kcol <= N_det_beta_unique) + + tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) + tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) + + if (kcol /= kcol_prev) then !if we've moved to a new unique beta determinant + call get_list_hp_banned_spin(tmp_det,N_hp,exc_is_banned_b1,spin_hp,sign_hp,idx_hp,2,$N_int,all_banned_b1) + if (all_banned_b1) then + kcol_prev = kcol + cycle + else ! get all unique beta dets connected to this one by a single excitation + call get_all_spin_singles_$N_int( & + psi_det_beta_unique, idx0, & + tmp_det(1,2), N_det_beta_unique, & + singles_b, n_singles_b) + kcol_prev = kcol + endif + else + if (all_banned_b1) cycle + endif + + ! at least some beta allowed + ! check alpha + call get_list_hp_banned_spin(tmp_det,N_hp,exc_is_banned_a1,spin_hp,sign_hp,idx_hp,1,$N_int,all_banned_a1) + if (all_banned_a1) cycle + + all_banned_ab1=.True. + do ii=1,N_hp + exc_is_banned_ab1(ii)=(exc_is_banned_a1(ii).or.exc_is_banned_b1(ii)) + all_banned_ab1 = (all_banned_ab1.and.exc_is_banned_ab1(ii)) + enddo + if (all_banned_ab1) cycle +! kcol_prev = kcol ! keep track of old col to see when we've moved to a new one + + ! Loop over singly excited beta columns + ! ------------------------------------- + + do i=1,n_singles_b ! loop over other columns in this row + lcol = singles_b(i) + + tmp_det2(1:$N_int,2) = psi_det_beta_unique(1:$N_int, lcol) + + call get_list_hp_banned_spin(tmp_det2,N_hp,exc_is_banned_b2,spin_hp,sign_hp,idx_hp,2,$N_int,all_banned_b2) + if (all_banned_b2) cycle + + l_a = psi_bilinear_matrix_columns_loc(lcol) ! location of start of this column within psi_bilinear_mat + ASSERT (l_a <= N_det) + + do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - l_a ! loop over rows in this column + lrow = psi_bilinear_matrix_rows(l_a) ! get row (index of unique alpha det) + ASSERT (lrow <= N_det_alpha_unique) + + buffer(1:$N_int,j) = psi_det_alpha_unique(1:$N_int, lrow) ! get alpha det + + ASSERT (l_a <= N_det) + idx(j) = l_a ! indices of dets within psi_bilinear_mat + l_a = l_a+1 + enddo + j = j-1 + ! get all alpha dets in this column that are connected to ref alpha by a single exc. + call get_all_spin_singles_$N_int( & + buffer, idx, tmp_det(1,1), j, & + singles_a, n_singles_a ) + + ! Loop over alpha singles + ! ----------------------- + + do k = 1,n_singles_a + l_a = singles_a(k) + ASSERT (l_a <= N_det) + + lrow = psi_bilinear_matrix_rows(l_a) + ASSERT (lrow <= N_det_alpha_unique) + + tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow) + call get_list_hp_banned_spin(tmp_det2,N_hp,exc_is_banned_a2,spin_hp,sign_hp,idx_hp,1,$N_int,all_banned_a2) + if (all_banned_a2) cycle + all_banned_ab12 = .True. + do ii=1,N_hp + exc_is_banned_ab12(ii)=((exc_is_banned_ab1(ii).or.exc_is_banned_b2(ii)).or.exc_is_banned_a2(ii)) + allowed_hp(ii)=(.not.exc_is_banned_ab12(ii)) + all_banned_ab12 = (all_banned_ab12.and.exc_is_banned_ab12(ii)) + enddo + if (all_banned_ab12) cycle + call i_h_j_double_alpha_beta_hp(tmp_det,tmp_det2,$N_int,hij_hp,N_hp,spin_hp,sign_hp,idx_hp,allowed_hp) + do l=1,N_hp + v_t(l,k_a) = v_t(l,k_a) + hij_hp(l) * u_t(l,l_a) + enddo + enddo + enddo + enddo + !$OMP END DO + + !$OMP DO SCHEDULE(dynamic,64) + do k_a=istart+ishift,iend,istep + + + ! Single and double alpha excitations + ! =================================== + + + ! Initial determinant is at k_a in alpha-major representation + ! ----------------------------------------------------------------------- + + krow = psi_bilinear_matrix_rows(k_a) + ASSERT (krow <= N_det_alpha_unique) + + kcol = psi_bilinear_matrix_columns(k_a) + ASSERT (kcol <= N_det_beta_unique) + + tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) + tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) + call get_list_hp_banned_ab(tmp_det,N_hp,exc_is_banned_ab1,spin_hp,sign_hp,idx_hp,$N_int,all_banned_ab1) + if (all_banned_ab1) cycle + + + ! Initial determinant is at k_b in beta-major representation + ! ---------------------------------------------------------------------- + + k_b = psi_bilinear_matrix_order_transp_reverse(k_a) + ASSERT (k_b <= N_det) + + spindet(1:$N_int) = tmp_det(1:$N_int,1) + + ! Loop inside the beta column to gather all the connected alphas + lcol = psi_bilinear_matrix_columns(k_a) + l_a = psi_bilinear_matrix_columns_loc(lcol) + do i=1,N_det_alpha_unique + if (l_a > N_det) exit + lcol = psi_bilinear_matrix_columns(l_a) + if (lcol /= kcol) exit + lrow = psi_bilinear_matrix_rows(l_a) + ASSERT (lrow <= N_det_alpha_unique) + + buffer(1:$N_int,i) = psi_det_alpha_unique(1:$N_int, lrow) + idx(i) = l_a + l_a = l_a+1 + enddo + i = i-1 + + !call get_all_spin_singles_and_doubles_$N_int( & + ! buffer, idx, spindet, i, & + ! singles_a, doubles, n_singles_a, n_doubles ) + call get_all_spin_singles_and_doubles( & + buffer, idx, spindet, $N_int, i, & + singles_a, doubles, n_singles_a, n_doubles ) + + ! Compute Hij for all alpha singles + ! ---------------------------------- + + tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) + do i=1,n_singles_a + l_a = singles_a(i) + ASSERT (l_a <= N_det) + + lrow = psi_bilinear_matrix_rows(l_a) + ASSERT (lrow <= N_det_alpha_unique) + + tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow) + call get_list_hp_banned_spin(tmp_det2,N_hp,exc_is_banned_a2,spin_hp,sign_hp,idx_hp,1,$N_int,all_banned_a2) + if (all_banned_a2) cycle + all_banned_ab12 = .True. + do ii=1,N_hp + exc_is_banned_ab12(ii)=(exc_is_banned_ab1(ii).or.exc_is_banned_a2(ii)) + allowed_hp(ii)=(.not.exc_is_banned_ab12(ii)) + all_banned_ab12 = (all_banned_ab12.and.exc_is_banned_ab12(ii)) + enddo + if (all_banned_ab12) cycle + call i_h_j_mono_spin_hp(tmp_det,tmp_det2,$N_int,1, hij_hp,N_hp,spin_hp,sign_hp,idx_hp,allowed_hp) + + do l=1,N_hp + v_t(l,k_a) = v_t(l,k_a) + hij_hp(l) * u_t(l,l_a) + ! single => sij = 0 + enddo + enddo + + ! Compute Hij for all alpha doubles + ! ---------------------------------- + + do i=1,n_doubles + l_a = doubles(i) + ASSERT (l_a <= N_det) + + lrow = psi_bilinear_matrix_rows(l_a) + ASSERT (lrow <= N_det_alpha_unique) + + call get_list_hp_banned_single_spin(psi_det_alpha_unique(1,lrow),N_hp,exc_is_banned_a2,spin_hp,sign_hp,idx_hp,1,$N_int,all_banned_a2) + if (all_banned_a2) cycle + all_banned_ab12 = .True. + do ii=1,N_hp + exc_is_banned_ab12(ii)=(exc_is_banned_ab1(ii).or.exc_is_banned_a2(ii)) + allowed_hp(ii)=(.not.exc_is_banned_ab12(ii)) + all_banned_ab12 = (all_banned_ab12.and.exc_is_banned_ab12(ii)) + enddo + if (all_banned_ab12) cycle + call i_h_j_double_spin_hp( tmp_det(1,1), psi_det_alpha_unique(1, lrow), $N_int,1,hij_hp,N_hp,spin_hp,sign_hp,idx_hp,allowed_hp) + do l=1,N_hp + v_t(l,k_a) = v_t(l,k_a) + hij_hp(l) * u_t(l,l_a) + ! same spin => sij = 0 + enddo + enddo + + + ! Single and double beta excitations + ! ================================== + + + ! Initial determinant is at k_a in alpha-major representation + ! ----------------------------------------------------------------------- + + krow = psi_bilinear_matrix_rows(k_a) + kcol = psi_bilinear_matrix_columns(k_a) + + tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) + tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) + + !! should already be done from top of loop? + !call get_list_hp_banned_ab(tmp_det,N_hp,exc_is_banned_ab1,spin_hp,sign_hp,idx_hp,$N_int,all_banned_ab1) + !if (all_banned_ab1) cycle + + spindet(1:$N_int) = tmp_det(1:$N_int,2) + + ! Initial determinant is at k_b in beta-major representation + ! ----------------------------------------------------------------------- + + k_b = psi_bilinear_matrix_order_transp_reverse(k_a) + ASSERT (k_b <= N_det) + + ! Loop inside the alpha row to gather all the connected betas + lrow = psi_bilinear_matrix_transp_rows(k_b) + l_b = psi_bilinear_matrix_transp_rows_loc(lrow) + do i=1,N_det_beta_unique + if (l_b > N_det) exit + lrow = psi_bilinear_matrix_transp_rows(l_b) + if (lrow /= krow) exit + lcol = psi_bilinear_matrix_transp_columns(l_b) + ASSERT (lcol <= N_det_beta_unique) + + buffer(1:$N_int,i) = psi_det_beta_unique(1:$N_int, lcol) + idx(i) = l_b + l_b = l_b+1 + enddo + i = i-1 + + !call get_all_spin_singles_and_doubles_$N_int( & + ! buffer, idx, spindet, i, & + ! singles_b, doubles, n_singles_b, n_doubles ) + call get_all_spin_singles_and_doubles( & + buffer, idx, spindet, $N_int, i, & + singles_b, doubles, n_singles_b, n_doubles ) + + ! Compute Hij for all beta singles + ! ---------------------------------- + + tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) + do i=1,n_singles_b + l_b = singles_b(i) + ASSERT (l_b <= N_det) + + lcol = psi_bilinear_matrix_transp_columns(l_b) + ASSERT (lcol <= N_det_beta_unique) + + tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, lcol) + call get_list_hp_banned_spin(tmp_det2,N_hp,exc_is_banned_b2,spin_hp,sign_hp,idx_hp,2,$N_int,all_banned_b2) + if (all_banned_b2) cycle + all_banned_ab12 = .True. + do ii=1,N_hp + exc_is_banned_ab12(ii)=(exc_is_banned_ab1(ii).or.exc_is_banned_b2(ii)) + allowed_hp(ii)=(.not.exc_is_banned_ab12(ii)) + all_banned_ab12 = (all_banned_ab12.and.exc_is_banned_ab12(ii)) + enddo + if (all_banned_ab12) cycle + call i_h_j_mono_spin_hp(tmp_det,tmp_det2,$N_int,2, hij_hp,N_hp,spin_hp,sign_hp,idx_hp,allowed_hp) + l_a = psi_bilinear_matrix_transp_order(l_b) + ASSERT (l_a <= N_det) + do l=1,N_hp + v_t(l,k_a) = v_t(l,k_a) + hij_hp(l) * u_t(l,l_a) + ! single => sij = 0 + enddo + enddo + + ! Compute Hij for all beta doubles + ! ---------------------------------- + + do i=1,n_doubles + l_b = doubles(i) + ASSERT (l_b <= N_det) + + lcol = psi_bilinear_matrix_transp_columns(l_b) + ASSERT (lcol <= N_det_beta_unique) + + call get_list_hp_banned_single_spin(psi_det_beta_unique(1,lcol),N_hp,exc_is_banned_b2,spin_hp,sign_hp,idx_hp,2,$N_int,all_banned_b2) + if (all_banned_b2) cycle + all_banned_ab12 = .True. + do ii=1,N_hp + exc_is_banned_ab12(ii)=(exc_is_banned_ab1(ii).or.exc_is_banned_b2(ii)) + allowed_hp(ii)=(.not.exc_is_banned_ab12(ii)) + all_banned_ab12 = (all_banned_ab12.and.exc_is_banned_ab12(ii)) + enddo + if (all_banned_ab12) cycle + call i_h_j_double_spin_hp( tmp_det(1,2), psi_det_beta_unique(1, lcol), $N_int,2,hij_hp,N_hp,spin_hp,sign_hp,idx_hp,allowed_hp) + l_a = psi_bilinear_matrix_transp_order(l_b) + ASSERT (l_a <= N_det) + + do l=1,N_hp + v_t(l,k_a) = v_t(l,k_a) + hij_hp(l) * u_t(l,l_a) + ! same spin => sij = 0 + enddo + enddo + + + ! Diagonal contribution + ! ===================== + + + ! Initial determinant is at k_a in alpha-major representation + ! ----------------------------------------------------------------------- + + krow = psi_bilinear_matrix_rows(k_a) + ASSERT (krow <= N_det_alpha_unique) + + kcol = psi_bilinear_matrix_columns(k_a) + ASSERT (kcol <= N_det_beta_unique) + + tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) + tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) + + call get_list_hp_banned_ab(tmp_det,N_hp,exc_is_banned_ab1,spin_hp,sign_hp,idx_hp,$N_int,all_banned_ab1) + if (all_banned_ab1) cycle + + double precision, external :: diag_H_mat_elem, diag_S_mat_elem + hii = diag_h_mat_elem(tmp_det,$N_int) + + do ii=1,N_hp + if(exc_is_banned_ab1(ii)) then + hii_hp(ii)=0.d0 + else + tmp_det2=tmp_det + na=elec_num_tab(spin_hp(ii)) + nb=elec_num_tab(iand(spin_hp(ii),1)+1) + hii_hp(ii)=hii + if (sign_hp(ii)>0) then + call ac_operator(idx_hp(ii),spin_hp(ii),tmp_det2,hii_hp(ii),$N_int,na,nb) + else + call a_operator(idx_hp(ii),spin_hp(ii),tmp_det2,hii_hp(ii),$N_int,na,nb) + endif + endif + v_t(ii,k_a) = v_t(ii,k_a) + (nuclear_repulsion + hii_hp(ii)) * u_t(ii,k_a) + enddo + + + end do + !$OMP END DO + deallocate(buffer, singles_a, singles_b, doubles, idx, & + exc_is_banned_a1, & + exc_is_banned_b1, & + exc_is_banned_a2, & + exc_is_banned_b2, & + exc_is_banned_ab1, & + exc_is_banned_ab12, & + allowed_hp, & + hij_hp, hii_hp ) + !$OMP END PARALLEL + deallocate(idx0) +end + +SUBST [ N_int ] + +1;; +2;; +3;; +4;; +N_int;; + +END_TEMPLATE + + + +subroutine i_h_j_double_spin_hp(key_i,key_j,Nint,ispin,hij_hp,N_hp,spin_hp,sign_hp,idx_hp,allowed_hp) + use bitmasks + implicit none + BEGIN_DOC + ! todo: maybe make new get_double_excitation_spin? + ! the 4 index ordering is already done in there, so we could avoid duplicating that work + ! Returns where i and j are determinants differing by a same-spin double excitation + END_DOC + integer, intent(in) :: Nint,ispin,N_hp + integer(bit_kind), intent(in) :: key_i(Nint), key_j(Nint) + complex*16, intent(out) :: hij_hp(N_hp) + integer, intent(in) :: spin_hp(N_hp), idx_hp(N_hp) + double precision, intent(in) :: sign_hp(N_hp) + logical, intent(in) :: allowed_hp(N_hp) + complex*16 :: hij0 + double precision :: phase_hp(N_hp) + integer :: exc(0:2,2) + double precision :: phase + complex*16, external :: get_mo_bielec_integral + integer :: i1,i2,i3,i4,j2,j3,ii + + PROVIDE big_array_exchange_integrals mo_bielec_integrals_in_map + + call get_double_excitation_spin(key_i,key_j,exc,phase,Nint) + hij0 = phase*(get_mo_bielec_integral( & + exc(1,1), & + exc(2,1), & + exc(1,2), & + exc(2,2), mo_integrals_map) - & + get_mo_bielec_integral( & + exc(1,1), & + exc(2,1), & + exc(2,2), & + exc(1,2), mo_integrals_map) ) + + ASSERT (exc(1,1) < exc(2,1)) + ASSERT (exc(1,2) < exc(2,2)) + i1=min(exc(1,1),exc(1,2)) + j2=max(exc(1,1),exc(1,2)) + j3=min(exc(2,1),exc(2,2)) + i4=max(exc(2,1),exc(2,2)) + i2=min(j2,j3) + i3=max(j2,j3) + + do ii=1,N_hp + if (allowed_hp(ii)) then + if (ispin.eq.spin_hp(ii)) then + if ((idx_hp(ii).lt.i1).or.(idx_hp(ii).gt.i4)) then + phase_hp(ii)=1.d0 + else if ((idx_hp(ii).lt.i2).or.(idx_hp(ii).gt.i3)) then + phase_hp(ii)=-1.d0 + else + phase_hp(ii)=1.d0 + endif + else + phase_hp(ii)=1.d0 + endif + else + phase_hp(ii)=0.d0 + endif + hij_hp(ii) = hij0 * phase_hp(ii) + enddo +end + +subroutine i_h_j_mono_spin_hp(key_i,key_j,Nint,spin,hij_hp,N_hp,spin_hp,sign_hp,idx_hp,allowed_hp) + use bitmasks + implicit none + BEGIN_DOC + ! todo: change this to use normal version of get_mono_excitation_from_fock + ! all info needed is in phase and hij, h/p part can happen after getting hij the normal way + ! Returns where i and j are determinants differing by a single excitation + END_DOC + integer, intent(in) :: Nint, spin, N_hp + integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) + complex*16, intent(out) :: hij_hp(N_hp) + integer, intent(in) :: spin_hp(N_hp), idx_hp(N_hp) + double precision, intent(in) :: sign_hp(N_hp) + logical, intent(in) :: allowed_hp(N_hp) + !double precision :: phase_hp(N_hp) + complex*16 :: hij0 + + integer :: exc(0:2,2) + double precision :: phase + + PROVIDE big_array_exchange_integrals mo_bielec_integrals_in_map + + call get_mono_excitation_spin(key_i(1,spin),key_j(1,spin),exc,phase,Nint) + + call get_mono_excitation_from_fock_hp(key_i,key_j,exc(1,1),exc(1,2),spin,phase,N_hp,hij_hp,spin_hp,sign_hp,idx_hp,allowed_hp) +end + +subroutine get_mono_excitation_from_fock_hp(det_1,det_2,h,p,spin,phase,N_hp,hij_hp,spin_hp,sign_hp,idx_hp,allowed_hp) + use bitmasks + implicit none + integer,intent(in) :: h,p,spin,N_hp + double precision, intent(in) :: phase + integer(bit_kind), intent(in) :: det_1(N_int,2), det_2(N_int,2) + complex*16, intent(out) :: hij_hp(N_hp) + integer, intent(in) :: spin_hp(N_hp), idx_hp(N_hp) + double precision, intent(in) :: sign_hp(N_hp) + logical, intent(in) :: allowed_hp(N_hp) + double precision :: phase_hp(N_hp) + complex*16 :: hij0 + integer :: low,high + + integer(bit_kind) :: differences(N_int,2) + integer(bit_kind) :: hole(N_int,2) + integer(bit_kind) :: partcl(N_int,2) + integer :: occ_hole(N_int*bit_kind_size,2) + integer :: occ_partcl(N_int*bit_kind_size,2) + integer :: n_occ_ab_hole(2),n_occ_ab_partcl(2) + integer :: i0,i,ii + do i = 1, N_int + differences(i,1) = xor(det_1(i,1),ref_closed_shell_bitmask(i,1)) + differences(i,2) = xor(det_1(i,2),ref_closed_shell_bitmask(i,2)) + hole(i,1) = iand(differences(i,1),ref_closed_shell_bitmask(i,1)) + hole(i,2) = iand(differences(i,2),ref_closed_shell_bitmask(i,2)) + partcl(i,1) = iand(differences(i,1),det_1(i,1)) + partcl(i,2) = iand(differences(i,2),det_1(i,2)) + enddo + call bitstring_to_list_ab(hole, occ_hole, n_occ_ab_hole, N_int) + call bitstring_to_list_ab(partcl, occ_partcl, n_occ_ab_partcl, N_int) + hij0 = fock_operator_closed_shell_ref_bitmask(h,p) + ! holes :: direct terms + do i0 = 1, n_occ_ab_hole(1) + i = occ_hole(i0,1) + hij0 -= big_array_coulomb_integrals(i,h,p) ! get_mo_bielec_integral_schwartz(h,i,p,i,mo_integrals_map) + enddo + do i0 = 1, n_occ_ab_hole(2) + i = occ_hole(i0,2) + hij0 -= big_array_coulomb_integrals(i,h,p) !get_mo_bielec_integral_schwartz(h,i,p,i,mo_integrals_map) + enddo + + ! holes :: exchange terms + do i0 = 1, n_occ_ab_hole(spin) + i = occ_hole(i0,spin) + hij0 += big_array_exchange_integrals(i,h,p) ! get_mo_bielec_integral_schwartz(h,i,i,p,mo_integrals_map) + enddo + + ! particles :: direct terms + do i0 = 1, n_occ_ab_partcl(1) + i = occ_partcl(i0,1) + hij0 += big_array_coulomb_integrals(i,h,p)!get_mo_bielec_integral_schwartz(h,i,p,i,mo_integrals_map) + enddo + do i0 = 1, n_occ_ab_partcl(2) + i = occ_partcl(i0,2) + hij0 += big_array_coulomb_integrals(i,h,p) !get_mo_bielec_integral_schwartz(h,i,p,i,mo_integrals_map) + enddo + + ! particles :: exchange terms + do i0 = 1, n_occ_ab_partcl(spin) + i = occ_partcl(i0,spin) + hij0 -= big_array_exchange_integrals(i,h,p)!get_mo_bielec_integral_schwartz(h,i,i,p,mo_integrals_map) + enddo + + low=min(h,p) + high=max(h,p) + +!! do ii=1,N_hp +!! if (.not.allowed_hp(ii)) then +!! phase_hp(ii) = 0.d0 +!! cycle +!! else if (spin_hp(ii).ne.spin) then +!! phase_hp(ii) = 1.d0 +!! else +!! if ((low.lt.idx_hp(ii)).and.(high.gt.idx_hp(ii))) then +!! phase_hp(ii) = -1.d0 +!! else +!! phase_hp(ii) = 1.d0 +!! endif +!! endif +!! enddo +!! +!! do ii=1,N_hp +!! if (allowed_hp(ii)) then +!! hij_hp(ii) = hij + sign_hp(ii) * big_array_coulomb_integrals(idx_hp(ii),h,p) +!! if (spin.eq.spin_hp(ii)) then +!! hij_hp(ii) = hij_hp(ii) - sign_hp(ii) * big_array_exchange_integrals(idx_hp(ii),h,p) +!! endif +!! else +!! hij_hp(ii) = 0.d0 +!! endif +!! enddo +!! +!! do ii=1,N_hp +!! hij_hp(ii) = hij_hp(ii) * phase_hp(ii) * phase +!! enddo + + do ii=1,N_hp + if (.not.allowed_hp(ii)) then + phase_hp(ii) = 0.d0 + hij_hp(ii) = 0.d0 + cycle + else if (spin.eq.spin_hp(ii)) then + hij_hp(ii) = hij0 + sign_hp(ii) *(big_array_coulomb_integrals(idx_hp(ii),h,p) - big_array_exchange_integrals(idx_hp(ii),h,p)) + if ((low.lt.idx_hp(ii)).and.(high.gt.idx_hp(ii))) then + phase_hp(ii) = -1.d0 + else + phase_hp(ii) = 1.d0 + endif + else + phase_hp(ii) = 1.d0 + hij_hp(ii) = hij0 + sign_hp(ii) * big_array_coulomb_integrals(idx_hp(ii),h,p) + endif + hij_hp(ii) = hij_hp(ii) * phase * phase_hp(ii) + enddo + +end + + +subroutine i_H_j_double_alpha_beta_hp(key_i,key_j,Nint,hij_hp,N_hp,spin_hp,sign_hp,idx_hp,allowed_hp) + use bitmasks + implicit none + BEGIN_DOC + ! Returns where i and j are determinants differing by an opposite-spin double excitation + END_DOC + integer, intent(in) :: Nint,N_hp + integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) + complex*16, intent(out) :: hij_hp(N_hp) + complex*16 :: hij0 + integer, intent(in) :: spin_hp(N_hp), idx_hp(N_hp) + double precision, intent(in) :: sign_hp(N_hp) + logical, intent(in) :: allowed_hp(N_hp) + double precision :: phase_hp(N_hp) + integer :: i + + integer :: lowhigh(2,2) + integer :: exc(0:2,2,2) + double precision :: phase, phase2 + complex*16, external :: get_mo_bielec_integral + + PROVIDE big_array_exchange_integrals mo_bielec_integrals_in_map + + call get_mono_excitation_spin(key_i(1,1),key_j(1,1),exc(0,1,1),phase,Nint) + call get_mono_excitation_spin(key_i(1,2),key_j(1,2),exc(0,1,2),phase2,Nint) + phase = phase*phase2 + + if (exc(1,1,1) == exc(1,2,2)) then + hij0 = big_array_exchange_integrals(exc(1,1,1),exc(1,1,2),exc(1,2,1)) + else if (exc(1,2,1) == exc(1,1,2)) then + hij0 = big_array_exchange_integrals(exc(1,2,1),exc(1,1,1),exc(1,2,2)) + else + hij0 = get_mo_bielec_integral( & + exc(1,1,1), & + exc(1,1,2), & + exc(1,2,1), & + exc(1,2,2) ,mo_integrals_map) + endif + + !todo: clean this up + ! if new particle/hole is between p/h of single exc of same spin, then parity changes, otherwise stays the same + ! value of Hij for double excitation is unchanged (new p/h is not one of the indices involved in the excitation) + + lowhigh(1,1)=min(exc(1,1,1),exc(1,2,1)) + lowhigh(2,1)=max(exc(1,1,1),exc(1,2,1)) + lowhigh(1,2)=min(exc(1,1,2),exc(1,2,2)) + lowhigh(2,2)=max(exc(1,1,2),exc(1,2,2)) + do i=1,N_hp + if (.not.allowed_hp(i)) then + phase_hp(i)=0.d0 + else if ((idx_hp(i).gt.lowhigh(1,spin_hp(i))).and.(idx_hp(i).lt.lowhigh(2,spin_hp(i)))) then + phase_hp(i)=-1.d0 + else + phase_hp(i)=1.d0 + endif + hij_hp(i)=hij0*phase*phase_hp(i) + enddo +end diff --git a/src/green/hu0_lanczos.irp.f b/src/green/hu0_lanczos.irp.f new file mode 100644 index 00000000..e4da5c78 --- /dev/null +++ b/src/green/hu0_lanczos.irp.f @@ -0,0 +1,405 @@ +! modified from H_S2_u_0_nstates_openmp in Davidson/u0Hu0.irp.f + +subroutine H_u_0_openmp(v_0,u_0,sze) + use bitmasks + implicit none + BEGIN_DOC + ! Computes v_0 = H|u_0> + ! + ! Assumes that the determinants are in psi_det + ! + ! istart, iend, ishift, istep are used in ZMQ parallelization. + END_DOC + integer :: N_st=1 + integer, intent(in) :: sze + complex*16, intent(inout) :: v_0(sze), u_0(sze) + integer :: k + call cdset_order(u_0(1),psi_bilinear_matrix_order,N_det) + v_0 = (0.d0,0.d0) + + call h_u_0_openmp_work(v_0,u_0,sze,1,N_det,0,1) + + call cdset_order(v_0(1),psi_bilinear_matrix_order_reverse,N_det) + call cdset_order(u_0(1),psi_bilinear_matrix_order_reverse,N_det) + +end + + +subroutine H_u_0_openmp_work(v_t,u_t,sze,istart,iend,ishift,istep) + use bitmasks + implicit none + BEGIN_DOC + ! Computes v_t = H|u_t> + ! + ! Default should be 1,N_det,0,1 + END_DOC + integer :: N_st=1 + integer, intent(in) :: sze,istart,iend,ishift,istep + complex*16, intent(in) :: u_t(N_det) + complex*16, intent(out) :: v_t(sze) + + + PROVIDE ref_bitmask_energy N_int + + select case (N_int) + case (1) + call H_u_0_openmp_work_1(v_t,u_t,sze,istart,iend,ishift,istep) + case (2) + call H_u_0_openmp_work_2(v_t,u_t,sze,istart,iend,ishift,istep) + case (3) + call H_u_0_openmp_work_3(v_t,u_t,sze,istart,iend,ishift,istep) + case (4) + call H_u_0_openmp_work_4(v_t,u_t,sze,istart,iend,ishift,istep) + case default + call H_u_0_openmp_work_N_int(v_t,u_t,sze,istart,iend,ishift,istep) + end select +end +BEGIN_TEMPLATE + +subroutine H_u_0_openmp_work_$N_int(v_t,u_t,sze,istart,iend,ishift,istep) + use bitmasks + implicit none + BEGIN_DOC + ! Computes v_t = H|u_t> + ! + ! Default should be 1,N_det,0,1 + END_DOC + integer :: N_st=1 + integer, intent(in) :: sze,istart,iend,ishift,istep + complex*16, intent(in) :: u_t(N_det) + complex*16, intent(out) :: v_t(sze) + + complex*16 :: hij + double precision :: hii + integer :: i,j,k,l + integer :: k_a, k_b, l_a, l_b, m_a, m_b + integer :: istate + integer :: krow, kcol, krow_b, kcol_b + integer :: lrow, lcol + integer :: mrow, mcol + integer(bit_kind) :: spindet($N_int) + integer(bit_kind) :: tmp_det($N_int,2) + integer(bit_kind) :: tmp_det2($N_int,2) + integer(bit_kind) :: tmp_det3($N_int,2) + integer(bit_kind), allocatable :: buffer(:,:) + integer :: n_doubles + integer, allocatable :: doubles(:) + integer, allocatable :: singles_a(:) + integer, allocatable :: singles_b(:) + integer, allocatable :: idx(:), idx0(:) + integer :: maxab, n_singles_a, n_singles_b, kcol_prev + integer*8 :: k8 + + maxab = max(N_det_alpha_unique, N_det_beta_unique)+1 + allocate(idx0(maxab)) + + do i=1,maxab + idx0(i) = i + enddo + + ! Prepare the array of all alpha single excitations + ! ------------------------------------------------- + + PROVIDE N_int nthreads_davidson + !$OMP PARALLEL DEFAULT(NONE) NUM_THREADS(nthreads_davidson) & + !$OMP SHARED(psi_bilinear_matrix_rows, N_det, & + !$OMP psi_bilinear_matrix_columns, & + !$OMP psi_det_alpha_unique, psi_det_beta_unique, & + !$OMP n_det_alpha_unique, n_det_beta_unique, N_int, & + !$OMP psi_bilinear_matrix_transp_rows, & + !$OMP psi_bilinear_matrix_transp_columns, & + !$OMP psi_bilinear_matrix_transp_order, N_st, & + !$OMP psi_bilinear_matrix_order_transp_reverse, & + !$OMP psi_bilinear_matrix_columns_loc, & + !$OMP psi_bilinear_matrix_transp_rows_loc, & + !$OMP istart, iend, istep, irp_here, v_t, & + !$OMP ishift, idx0, u_t, maxab) & + !$OMP PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i, & + !$OMP lcol, lrow, l_a, l_b, & + !$OMP buffer, doubles, n_doubles, & + !$OMP tmp_det2, hii, hij, idx, l, kcol_prev, & + !$OMP singles_a, n_singles_a, singles_b, & + !$OMP n_singles_b, k8) + + ! Alpha/Beta double excitations + ! ============================= + + allocate( buffer($N_int,maxab), & + singles_a(maxab), & + singles_b(maxab), & + doubles(maxab), & + idx(maxab)) + + kcol_prev=-1 + + ASSERT (iend <= N_det) + ASSERT (istart > 0) + ASSERT (istep > 0) + + !$OMP DO SCHEDULE(dynamic,64) + do k_a=istart+ishift,iend,istep + + krow = psi_bilinear_matrix_rows(k_a) + ASSERT (krow <= N_det_alpha_unique) + + kcol = psi_bilinear_matrix_columns(k_a) + ASSERT (kcol <= N_det_beta_unique) + + tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) + tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) + + if (kcol /= kcol_prev) then + call get_all_spin_singles_$N_int( & + psi_det_beta_unique, idx0, & + tmp_det(1,2), N_det_beta_unique, & + singles_b, n_singles_b) + endif + kcol_prev = kcol + + ! Loop over singly excited beta columns + ! ------------------------------------- + + do i=1,n_singles_b + lcol = singles_b(i) + + tmp_det2(1:$N_int,2) = psi_det_beta_unique(1:$N_int, lcol) + + l_a = psi_bilinear_matrix_columns_loc(lcol) + ASSERT (l_a <= N_det) + + do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - l_a + lrow = psi_bilinear_matrix_rows(l_a) + ASSERT (lrow <= N_det_alpha_unique) + + buffer(1:$N_int,j) = psi_det_alpha_unique(1:$N_int, lrow) + + ASSERT (l_a <= N_det) + idx(j) = l_a + l_a = l_a+1 + enddo + j = j-1 + + call get_all_spin_singles_$N_int( & + buffer, idx, tmp_det(1,1), j, & + singles_a, n_singles_a ) + + ! Loop over alpha singles + ! ----------------------- + + do k = 1,n_singles_a + l_a = singles_a(k) + ASSERT (l_a <= N_det) + + lrow = psi_bilinear_matrix_rows(l_a) + ASSERT (lrow <= N_det_alpha_unique) + + tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow) + call i_h_j_double_alpha_beta(tmp_det,tmp_det2,$N_int,hij) + v_t(k_a) = v_t(k_a) + hij * u_t(l_a) + enddo + enddo + enddo + !$OMP END DO + + !$OMP DO SCHEDULE(dynamic,64) + do k_a=istart+ishift,iend,istep + + + ! Single and double alpha excitations + ! =================================== + + + ! Initial determinant is at k_a in alpha-major representation + ! ----------------------------------------------------------------------- + + krow = psi_bilinear_matrix_rows(k_a) + ASSERT (krow <= N_det_alpha_unique) + + kcol = psi_bilinear_matrix_columns(k_a) + ASSERT (kcol <= N_det_beta_unique) + + tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) + tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) + + ! Initial determinant is at k_b in beta-major representation + ! ---------------------------------------------------------------------- + + k_b = psi_bilinear_matrix_order_transp_reverse(k_a) + ASSERT (k_b <= N_det) + + spindet(1:$N_int) = tmp_det(1:$N_int,1) + + ! Loop inside the beta column to gather all the connected alphas + lcol = psi_bilinear_matrix_columns(k_a) + l_a = psi_bilinear_matrix_columns_loc(lcol) + do i=1,N_det_alpha_unique + if (l_a > N_det) exit + lcol = psi_bilinear_matrix_columns(l_a) + if (lcol /= kcol) exit + lrow = psi_bilinear_matrix_rows(l_a) + ASSERT (lrow <= N_det_alpha_unique) + + buffer(1:$N_int,i) = psi_det_alpha_unique(1:$N_int, lrow) + idx(i) = l_a + l_a = l_a+1 + enddo + i = i-1 + + !call get_all_spin_singles_and_doubles_$N_int( & + ! buffer, idx, spindet, i, & + ! singles_a, doubles, n_singles_a, n_doubles ) + call get_all_spin_singles_and_doubles( & + buffer, idx, spindet, $N_int, i, & + singles_a, doubles, n_singles_a, n_doubles ) + + ! Compute Hij for all alpha singles + ! ---------------------------------- + + tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) + do i=1,n_singles_a + l_a = singles_a(i) + ASSERT (l_a <= N_det) + + lrow = psi_bilinear_matrix_rows(l_a) + ASSERT (lrow <= N_det_alpha_unique) + + tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow) + call i_H_j_mono_spin( tmp_det, tmp_det2, $N_int, 1, hij) + + v_t(k_a) = v_t(k_a) + hij * u_t(l_a) + enddo + + + ! Compute Hij for all alpha doubles + ! ---------------------------------- + + do i=1,n_doubles + l_a = doubles(i) + ASSERT (l_a <= N_det) + + lrow = psi_bilinear_matrix_rows(l_a) + ASSERT (lrow <= N_det_alpha_unique) + + call i_H_j_double_spin( tmp_det(1,1), psi_det_alpha_unique(1, lrow), $N_int, hij) + v_t(k_a) = v_t(k_a) + hij * u_t(l_a) + enddo + + + ! Single and double beta excitations + ! ================================== + + + ! Initial determinant is at k_a in alpha-major representation + ! ----------------------------------------------------------------------- + + krow = psi_bilinear_matrix_rows(k_a) + kcol = psi_bilinear_matrix_columns(k_a) + + tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) + tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) + + spindet(1:$N_int) = tmp_det(1:$N_int,2) + + ! Initial determinant is at k_b in beta-major representation + ! ----------------------------------------------------------------------- + + k_b = psi_bilinear_matrix_order_transp_reverse(k_a) + ASSERT (k_b <= N_det) + + ! Loop inside the alpha row to gather all the connected betas + lrow = psi_bilinear_matrix_transp_rows(k_b) + l_b = psi_bilinear_matrix_transp_rows_loc(lrow) + do i=1,N_det_beta_unique + if (l_b > N_det) exit + lrow = psi_bilinear_matrix_transp_rows(l_b) + if (lrow /= krow) exit + lcol = psi_bilinear_matrix_transp_columns(l_b) + ASSERT (lcol <= N_det_beta_unique) + + buffer(1:$N_int,i) = psi_det_beta_unique(1:$N_int, lcol) + idx(i) = l_b + l_b = l_b+1 + enddo + i = i-1 + + !call get_all_spin_singles_and_doubles_$N_int( & + ! buffer, idx, spindet, i, & + ! singles_b, doubles, n_singles_b, n_doubles ) + call get_all_spin_singles_and_doubles( & + buffer, idx, spindet, $N_int, i, & + singles_b, doubles, n_singles_b, n_doubles ) + + ! Compute Hij for all beta singles + ! ---------------------------------- + + tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) + do i=1,n_singles_b + l_b = singles_b(i) + ASSERT (l_b <= N_det) + + lcol = psi_bilinear_matrix_transp_columns(l_b) + ASSERT (lcol <= N_det_beta_unique) + + tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, lcol) + call i_H_j_mono_spin( tmp_det, tmp_det2, $N_int, 2, hij) + l_a = psi_bilinear_matrix_transp_order(l_b) + ASSERT (l_a <= N_det) + v_t(k_a) = v_t(k_a) + hij * u_t(l_a) + enddo + + ! Compute Hij for all beta doubles + ! ---------------------------------- + + do i=1,n_doubles + l_b = doubles(i) + ASSERT (l_b <= N_det) + + lcol = psi_bilinear_matrix_transp_columns(l_b) + ASSERT (lcol <= N_det_beta_unique) + + call i_H_j_double_spin( tmp_det(1,2), psi_det_beta_unique(1, lcol), $N_int, hij) + l_a = psi_bilinear_matrix_transp_order(l_b) + ASSERT (l_a <= N_det) + + v_t(k_a) = v_t(k_a) + hij * u_t(l_a) + enddo + + + ! Diagonal contribution + ! ===================== + + + ! Initial determinant is at k_a in alpha-major representation + ! ----------------------------------------------------------------------- + + krow = psi_bilinear_matrix_rows(k_a) + ASSERT (krow <= N_det_alpha_unique) + + kcol = psi_bilinear_matrix_columns(k_a) + ASSERT (kcol <= N_det_beta_unique) + + tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) + tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) + + !double precision, external :: diag_H_mat_elem, diag_S_mat_elem + double precision, external :: diag_H_mat_elem + hii = diag_H_mat_elem(tmp_det,$N_int) + v_t(k_a) = v_t(k_a) + hii * u_t(k_a) + + end do + !$OMP END DO + deallocate(buffer, singles_a, singles_b, doubles, idx) + !$OMP END PARALLEL + +end + +SUBST [ N_int ] + +1;; +2;; +3;; +4;; +N_int;; + +END_TEMPLATE + diff --git a/src/green/lanczos.irp.f b/src/green/lanczos.irp.f new file mode 100644 index 00000000..a2557abb --- /dev/null +++ b/src/green/lanczos.irp.f @@ -0,0 +1,882 @@ + + +BEGIN_PROVIDER [ integer, n_green_vec ] + implicit none + BEGIN_DOC + ! number of particles/holes to use for spectral density calc. + ! just set to 2 for now (homo and lumo) + END_DOC + n_green_vec = 2 +END_PROVIDER + + BEGIN_PROVIDER [ integer, green_idx, (n_green_vec) ] +&BEGIN_PROVIDER [ integer, green_idx_int, (n_green_vec) ] +&BEGIN_PROVIDER [ integer, green_idx_bit, (n_green_vec) ] +&BEGIN_PROVIDER [ integer, green_spin, (n_green_vec) ] +&BEGIN_PROVIDER [ double precision, green_sign, (n_green_vec) ] + implicit none + BEGIN_DOC + ! description of particles/holes to be used in spectral density calculation + ! green_idx: orbital index of particle/hole + ! green_idx_{int,bit}: location of idx within determinant bitstring + ! green_spin: 1(alpha) or 2(beta) + ! green_sign: 1(particle) or -1(hole) + END_DOC + integer :: s1,s2,i1,i2 + integer :: i + + integer :: idx_homo_lumo(2), spin_homo_lumo(2) + logical :: has_idx,has_spin,has_sign,has_lanc + integer :: nlanc + ! needs psi_det, mo_tot_num, N_int, mo_bielec_integral_jj, mo_mono_elec_integral_diag + call ezfio_has_green_green_idx(has_idx) + call ezfio_has_green_green_spin(has_spin) + call ezfio_has_green_green_sign(has_sign) +! call ezfio_has_green_n_lanczos_complete(has_lanc) + call ezfio_get_green_n_lanczos_complete(nlanc) + if (has_idx.and.has_spin.and.has_sign) then + print*,'reading idx,spin,sign' + call ezfio_get_green_green_idx(green_idx) + call ezfio_get_green_green_spin(green_spin) + call ezfio_get_green_green_sign(green_sign) + else if (nlanc.gt.0) then + stop 'problem with lanczos restart; need idx, spin, sign' + else + print*,'new lanczos calculation, finding homo/lumo' + call get_homo_lumo(psi_det(1:N_int,1:2,1),N_int,mo_tot_num,idx_homo_lumo,spin_homo_lumo) + + ! homo + green_idx(1)=idx_homo_lumo(1) + green_spin(1)=spin_homo_lumo(1) + green_sign(1)=-1.d0 + + ! lumo + green_idx(2)=idx_homo_lumo(2) + green_spin(2)=spin_homo_lumo(2) + green_sign(2)=1.d0 + + call ezfio_set_green_green_idx(green_idx) + call ezfio_set_green_green_spin(green_spin) + call ezfio_set_green_green_sign(green_sign) + endif + + + +! if (nlanc.gt.0) then +! ! call ezfio_get_green_n_lanczos_complete(nlanc) +! print*,'restarting from previous lanczos',nlanc +! if (has_idx.and.has_spin.and.has_sign) then +! print*,'reading idx,spin,sign' +! call ezfio_get_green_green_idx(green_idx) +! call ezfio_get_green_green_spin(green_spin) +! call ezfio_get_green_green_sign(green_sign) +! else +! stop 'problem with lanczos restart; need idx, spin, sign' +! endif +! else +! print*,'new lanczos calculation, finding homo/lumo' +! call get_homo_lumo(psi_det(1:N_int,1:2,1),N_int,mo_tot_num,idx_homo_lumo,spin_homo_lumo) +! +! ! homo +! green_idx(1)=idx_homo_lumo(1) +! green_spin(1)=spin_homo_lumo(1) +! green_sign(1)=-1.d0 +! +! ! lumo +! green_idx(2)=idx_homo_lumo(2) +! green_spin(2)=spin_homo_lumo(2) +! green_sign(2)=1.d0 +! +! call ezfio_set_green_green_idx(green_idx) +! call ezfio_set_green_green_spin(green_spin) +! call ezfio_set_green_green_sign(green_sign) +! endif + + do i=1,n_green_vec + call get_orb_int_bit(green_idx(i),green_idx_int(i),green_idx_bit(i)) + print*,i,green_idx(i),green_idx_int(i),green_idx_bit(i),green_spin(i),green_sign(i) + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, green_det_phase, (N_det,n_green_vec) ] + implicit none + BEGIN_DOC + ! for each det in psi, compute phase for each particle/hole excitation + ! each element should be +/-1 or 0 + END_DOC + integer :: i + double precision :: phase_tmp(n_green_vec) + PROVIDE psi_det green_idx + + do i=1,N_det + call get_phase_hp(green_idx_int,green_idx_bit,green_spin,green_sign,psi_det(1,1,i),phase_tmp,N_int,n_green_vec) + green_det_phase(i,1:n_green_vec) = phase_tmp(1:n_green_vec) + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ complex*16, u1_lanczos, (N_det,n_green_vec) ] + implicit none + BEGIN_DOC + ! initial lanczos vectors + ! must be normalized + END_DOC + + integer :: i,j + + do j=1,n_green_vec + do i=1,N_det + u1_lanczos(i,j)=green_det_phase(i,j)*psi_coef(i,1) + enddo + call normalize_complex(u1_lanczos(:,j),N_det) + enddo + +END_PROVIDER + +! BEGIN_PROVIDER [ double precision, alpha_lanczos, (n_green_vec,n_lanczos_iter) ] +!&BEGIN_PROVIDER [ double precision, beta_lanczos, (n_green_vec,n_lanczos_iter) ] + BEGIN_PROVIDER [ double precision, alpha_lanczos, (n_lanczos_iter,n_green_vec) ] +&BEGIN_PROVIDER [ double precision, beta_lanczos, (n_lanczos_iter,n_green_vec) ] +&BEGIN_PROVIDER [ complex*16, un_lanczos, (N_det,n_green_vec) ] +&BEGIN_PROVIDER [ complex*16, vn_lanczos, (N_det,n_green_vec) ] +&BEGIN_PROVIDER [ double precision, lanczos_eigvals, (n_lanczos_iter,n_green_vec) ] + implicit none + BEGIN_DOC + ! for each particle/hole: + ! provide alpha and beta for tridiagonal form of H + ! un, vn lanczos vectors from latest iteration + ! lanczos_eigvals: eigenvalues of tridiagonal form of H + END_DOC + PROVIDE lanczos_debug_print n_lanczos_debug + complex*16, allocatable :: work(:,:) +! double precision :: alpha_tmp,beta_tmp + double precision, allocatable :: alpha_tmp(:),beta_tmp(:) + double precision, allocatable :: alpha_tmp_vec(:,:), beta_tmp_vec(:,:) + integer :: i,j + integer :: n_lanc_new_tmp, n_lanc_old_tmp + call ezfio_get_green_n_lanczos_iter(n_lanc_new_tmp) + call ezfio_get_green_n_lanczos_complete(n_lanc_old_tmp) + + if ((n_lanczos_complete).gt.0) then +! allocate(alpha_tmp_vec(n_green_vec,n_lanczos_complete),beta_tmp_vec(n_green_vec,n_lanczos_complete)) + allocate(alpha_tmp_vec(n_lanczos_complete,n_green_vec),beta_tmp_vec(n_lanczos_complete,n_green_vec)) + logical :: has_un_lanczos, has_vn_lanczos + call ezfio_has_green_un_lanczos(has_un_lanczos) + call ezfio_has_green_vn_lanczos(has_vn_lanczos) + if (has_un_lanczos.and.has_vn_lanczos) then + call ezfio_get_green_un_lanczos(un_lanczos) + call ezfio_get_green_vn_lanczos(vn_lanczos) +! if (lanczos_debug_print) then +! print*,'uu,vv read from disk' +! do i=1,n_lanczos_debug +! write(6,'(4(E25.15))')un_lanczos(i),vn_lanczos(i) +! enddo +! endif + else + print*,'problem reading lanczos vectors for restart' + stop + endif + logical :: has_alpha_lanczos, has_beta_lanczos + call ezfio_has_green_alpha_lanczos(has_alpha_lanczos) + call ezfio_has_green_beta_lanczos(has_beta_lanczos) + if (has_alpha_lanczos.and.has_beta_lanczos) then + call ezfio_set_green_n_lanczos_iter(n_lanc_old_tmp) + call ezfio_get_green_alpha_lanczos(alpha_tmp_vec) + call ezfio_get_green_beta_lanczos(beta_tmp_vec) + call ezfio_set_green_n_lanczos_iter(n_lanc_new_tmp) + do j=1,n_green_vec + do i=1,n_lanczos_complete + alpha_lanczos(i,j)=alpha_tmp_vec(i,j) + beta_lanczos(i,j)=beta_tmp_vec(i,j) + enddo + enddo + else + print*,'problem reading lanczos alpha, beta for restart' + stop + endif + deallocate(alpha_tmp_vec,beta_tmp_vec) + else + call write_time(6) + print*,'no saved lanczos vectors. starting lanczos' + PROVIDE u1_lanczos + un_lanczos=u1_lanczos + allocate(work(N_det,n_green_vec),alpha_tmp(n_green_vec),beta_tmp(n_green_vec)) + call lanczos_h_init_hp(un_lanczos,vn_lanczos,work,N_det,alpha_tmp,beta_tmp,& + n_green_vec,green_spin,green_sign,green_idx) + do i=1,n_green_vec + alpha_lanczos(1,i)=alpha_tmp(i) + beta_lanczos(1,i)=beta_tmp(i) + enddo + n_lanczos_complete=1 + deallocate(work,alpha_tmp,beta_tmp) + endif + + allocate(work(N_det,n_green_vec),alpha_tmp(n_green_vec),beta_tmp(n_green_vec)) + do i=n_lanczos_complete+1,n_lanczos_iter + call write_time(6) + print*,'starting lanczos iteration',i + call lanczos_h_step_hp(un_lanczos,vn_lanczos,work,N_det,alpha_tmp,beta_tmp,& + n_green_vec,green_spin,green_sign,green_idx) + do j=1,n_green_vec + alpha_lanczos(i,j)=alpha_tmp(j) + beta_lanczos(i,j)=beta_tmp(j) + enddo + n_lanczos_complete=n_lanczos_complete+1 + enddo + deallocate(work,alpha_tmp,beta_tmp) + + call ezfio_set_green_alpha_lanczos(alpha_lanczos) + call ezfio_set_green_beta_lanczos(beta_lanczos) + call ezfio_set_green_un_lanczos(un_lanczos) + call ezfio_set_green_vn_lanczos(vn_lanczos) + call ezfio_set_green_n_lanczos_complete(n_lanczos_complete) + + call diag_lanczos_vals_hp(alpha_lanczos, beta_lanczos, n_lanczos_complete, lanczos_eigvals,& + n_lanczos_iter,n_green_vec) + call ezfio_set_green_lanczos_eigvals(lanczos_eigvals) + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, delta_omega ] + implicit none + BEGIN_DOC + ! step size between frequency points for spectral density calculation + ! calculated from min, max, and number of steps + END_DOC + delta_omega=(omega_max-omega_min)/n_omega +END_PROVIDER + +BEGIN_PROVIDER [ double precision, omega_list, (n_omega) ] + implicit none + BEGIN_DOC + ! list of frequencies at which to compute spectral density + END_DOC + + integer :: i + double precision :: omega_i + PROVIDE delta_omega + do i=1,n_omega + omega_list(i) = omega_min + (i-1)*delta_omega + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, spectral_lanczos, (n_omega,n_green_vec) ] + implicit none + BEGIN_DOC + ! spectral density A(omega) calculated from lanczos alpha/beta + ! calculated for n_omega points between omega_min and omega_max + END_DOC + + integer :: i,j + double precision :: omega_i + complex*16 :: z_i + !double precision :: spec_lanc_rev + double precision :: spec_lanc_rev_sign + logical :: has_ci_energy + double precision :: ref_energy_0 + PROVIDE delta_omega alpha_lanczos beta_lanczos omega_list + call ezfio_has_full_ci_zmq_energy(has_ci_energy) + if (has_ci_energy) then + call ezfio_get_full_ci_zmq_energy(ref_energy_0) + else + print*,'no reference energy from full_ci_zmq, exiting' + stop + endif + + + do i=1,n_omega + omega_i = omega_list(i) + z_i = dcmplx(omega_i,gf_epsilon) + do j=1,n_green_vec +! spectral_lanczos(i,j) = spec_lanc_rev(n_lanczos_iter,alpha_lanczos(:,j),beta_lanczos(:,j),z_i) + spectral_lanczos(i,j) = spec_lanc_rev_sign(n_lanczos_iter, & + alpha_lanczos(:,j), & + beta_lanczos(:,j), & + z_i - green_sign(j)*ref_energy_0, & + green_sign(j)) + enddo + enddo + +END_PROVIDER + +double precision function spec_lanc(n_lanc_iter,alpha,beta,z) + include 'constants.include.F' + implicit none + BEGIN_DOC + ! input: + ! alpha, beta: from tridiagonal form of H (obtain via lanczos) + ! beta and alpha same size (beta(1) is not used) + ! n_lanc_iter: size of alpha, beta + ! z: omega + i*epsilon + ! omega is frequency for which spectral density is to be computed + ! epsilon is magnitude of infinitesimal imaginary term + ! output: + ! spec_lanc: spectral density A(omega) + ! + ! uses inv_pi=(1.d0/pi) from constants + END_DOC + integer, intent(in) :: n_lanc_iter + double precision, intent(in) :: alpha(n_lanc_iter), beta(n_lanc_iter) + complex*16, intent(in) :: z + + complex*16 bigAj2,bigAj1,bigAj0 + complex*16 bigBj2,bigBj1,bigBj0 + integer :: j + ! init for j=1 + ! bigAj2 is A(j-2) + ! bigAj1 is A(j-1) + ! etc. + + bigAj2=1.d0 ! A(-1) + bigAj1=0.d0 ! A(0) + bigAj0=1.d0 ! A(1) + + bigBj2=0.d0 ! B(-1) + bigBj1=1.d0 ! B(0) + bigBj0=z-alpha(1) ! B(1) + + do j=2,n_lanc_iter + bigAj2=bigAj1 + bigAj1=bigAj0 + bigAj0=(z-alpha(j))*bigAj1 - beta(j)**2*bigAj2 + + bigBj2=bigBj1 + bigBj1=bigBj0 + bigBj0=(z-alpha(j))*bigBj1 - beta(j)**2*bigBj2 + enddo + spec_lanc=-imag(bigAj0/bigBj0)*inv_pi +end + +double precision function spec_lanc_rev(n_lanc_iter,alpha,beta,z) + include 'constants.include.F' + implicit none + BEGIN_DOC + ! reverse iteration is more numerically stable + ! input: + ! alpha, beta: from tridiagonal form of H (obtain via lanczos) + ! beta and alpha same size (beta(1) is not used) + ! n_lanc_iter: size of alpha, beta + ! z: omega + i*epsilon + ! omega is frequency for which spectral density is to be computed + ! epsilon is magnitude of infinitesimal imaginary term + ! output: + ! spec_lanc: spectral density A(omega) + ! + ! uses inv_pi=(1.d0/pi) from constants + END_DOC + integer, intent(in) :: n_lanc_iter + double precision, intent(in) :: alpha(n_lanc_iter), beta(n_lanc_iter) + complex*16, intent(in) :: z + + complex*16 :: tmp + integer :: j + + tmp=(0.d0,0.d0) + do j=n_lanc_iter,2,-1 + tmp=-beta(j)**2/(z-alpha(j)+tmp) + enddo + tmp=1.d0/(z-alpha(1)+tmp) + spec_lanc_rev=-imag(tmp)*inv_pi +end + +double precision function spec_lanc_rev_sign(n_lanc_iter,alpha,beta,z,g_sign) + include 'constants.include.F' + implicit none + BEGIN_DOC + ! reverse iteration is more numerically stable + ! input: + ! alpha, beta: from tridiagonal form of H (obtain via lanczos) + ! beta and alpha same size (beta(1) is not used) + ! n_lanc_iter: size of alpha, beta + ! z: omega + i*epsilon + ! omega is frequency for which spectral density is to be computed + ! epsilon is magnitude of infinitesimal imaginary term + ! output: + ! spec_lanc: spectral density A(omega) + ! + ! uses inv_pi=(1.d0/pi) from constants + END_DOC + integer, intent(in) :: n_lanc_iter + double precision, intent(in) :: alpha(n_lanc_iter), beta(n_lanc_iter) + complex*16, intent(in) :: z + double precision, intent(in) :: g_sign + + complex*16 :: tmp + integer :: j + + tmp=(0.d0,0.d0) + do j=n_lanc_iter,2,-1 + tmp=-beta(j)**2/(z+g_sign*alpha(j)+tmp) + enddo + tmp=1.d0/(z+g_sign*alpha(1)+tmp) + spec_lanc_rev_sign=-imag(tmp)*inv_pi +end + + +subroutine lanczos_h_init_hp(uu,vv,work,sze,alpha_i,beta_i,ng,spin_hp,sign_hp,idx_hp) + implicit none + integer, intent(in) :: sze,ng + complex*16, intent(in) :: uu(sze,ng) + complex*16, intent(out) :: vv(sze,ng) + complex*16 :: work(sze,ng) + double precision, intent(out) :: alpha_i(ng), beta_i(ng) + integer, intent(in) :: spin_hp(ng), idx_hp(ng) + double precision, intent(in) :: sign_hp(ng) + + double precision, external :: dznrm2 + complex*16, external :: u_dot_v_complex + integer :: i,j + + BEGIN_DOC + ! initial step for lanczos tridiagonalization of H for multiple holes/particles + ! uu is array of initial vectors u1 (creation/annihilation operator applied to psi) + ! output vv is array of lanczos v1 (one for each hole/particle) + END_DOC + + print *,'starting lanczos' + print *,'sze = ',sze + + ! |uu> is |u(1)> + + ! |w(1)> = H|u(1)> + ! |work> is now |w(1)> + call compute_hu_hp(uu,work,ng,sze,spin_hp,sign_hp,idx_hp) + + ! alpha(n+1) = + do i=1,ng + alpha_i(i)=real(u_dot_v_complex(uu(1:sze,i),work(1:sze,i),sze)) + enddo + + do j=1,ng + do i=1,sze + vv(i,j)=work(i,j)-alpha_i(j)*uu(i,j) +! write(6,'(7(E25.15))')uu(i,j),vv(i,j),work(i,j),alpha_i(j) + enddo + enddo + + beta_i=0.d0 + ! |vv> is |v(1)> + ! |uu> is |u(1)> +end + +subroutine lanczos_h_step_hp(uu,vv,work,sze,alpha_i,beta_i,ng,spin_hp,sign_hp,idx_hp) + implicit none + integer, intent(in) :: sze,ng + complex*16, intent(inout) :: uu(sze,ng),vv(sze,ng) + complex*16, intent(out) :: work(sze,ng) + double precision, intent(out) :: alpha_i(ng), beta_i(ng) + integer, intent(in) :: spin_hp(ng), sign_hp(ng), idx_hp(ng) + + double precision, external :: dznrm2 + complex*16, external :: u_dot_v_complex + integer :: i,j + complex*16 :: tmp_c16 + BEGIN_DOC + ! lanczos tridiagonalization of H + ! n_lanc_iter is number of lanczos iterations + ! u1 is initial lanczos vector + ! u1 should be normalized + END_DOC + + ! |vv> is |v(n)> + ! |uu> is |u(n)> + + ! compute beta(n+1) + do j=1,ng + beta_i(j)=dznrm2(sze,vv(:,j),1) + ! |vv> is now |u(n+1)> + call zdscal(sze,(1.d0/beta_i(j)),vv(:,j),1) + enddo + + ! |w(n+1)> = H|u(n+1)> + ! |work> is now |w(n+1)> + call compute_hu_hp(vv,work,ng,sze,spin_hp,sign_hp,idx_hp) + + ! alpha(n+1) = + do i=1,ng + alpha_i(i)=real(u_dot_v_complex(vv(1:sze,i),work(1:sze,i),sze)) + enddo + + do j=1,ng + do i=1,sze + tmp_c16=work(i,j)-alpha_i(j)*vv(i,j)-beta_i(j)*uu(i,j) + uu(i,j)=vv(i,j) + vv(i,j)=tmp_c16 + enddo + enddo + ! |vv> is |v(n+1)> + ! |uu> is |u(n+1)> +end + + +subroutine lanczos_h_init(uu,vv,work,sze,alpha_i,beta_i) + implicit none + integer, intent(in) :: sze + complex*16, intent(inout) :: uu(sze) + complex*16, intent(out) :: vv(sze) + complex*16 :: work(sze) + double precision, intent(out) :: alpha_i, beta_i + + double precision, external :: dznrm2 + complex*16, external :: u_dot_v_complex + integer :: i + + BEGIN_DOC + ! lanczos tridiagonalization of H + ! n_lanc_iter is number of lanczos iterations + ! u1 is initial lanczos vector + ! u1 should be normalized + END_DOC + + print *,'starting lanczos' + print *,'sze = ',sze + ! exit if u1 is not normalized +! beta_norm = dznrm2(h_size,u1,1) +! if (dabs(beta_norm-1.d0) .gt. 1.d-6) then +! print *, 'Error: initial Lanczos vector is not normalized' +! stop -1 +! endif + + ! |uu> is |u(1)> + + ! |w(1)> = H|u(1)> + ! |work> is now |w(1)> + call compute_hu(uu,work,sze) + + ! alpha(n+1) = + alpha_i=real(u_dot_v_complex(uu,work,sze)) + + do i=1,sze + vv(i)=work(i)-alpha_i*uu(i) + enddo + beta_i=0.d0 + if (lanczos_debug_print) then + print*,'init uu,vv,work' + do i=1,n_lanczos_debug + write(6,'(6(E25.15))')uu(i),vv(i),work(i) + enddo + endif + ! |vv> is |v(1)> + ! |uu> is |u(1)> +end + +subroutine lanczos_h_step(uu,vv,work,sze,alpha_i,beta_i) + implicit none + integer, intent(in) :: sze + complex*16, intent(inout) :: uu(sze),vv(sze) + complex*16, intent(out) :: work(sze) + double precision, intent(out) :: alpha_i, beta_i + + double precision, external :: dznrm2 + complex*16, external :: u_dot_v_complex + integer :: i + complex*16 :: tmp_c16 + BEGIN_DOC + ! lanczos tridiagonalization of H + ! n_lanc_iter is number of lanczos iterations + ! u1 is initial lanczos vector + ! u1 should be normalized + END_DOC + + ! exit if u1 is not normalized +! beta_norm = dznrm2(h_size,u1,1) +! if (dabs(beta_norm-1.d0) .gt. 1.d-6) then +! print *, 'Error: initial Lanczos vector is not normalized' +! stop -1 +! endif + + ! |vv> is |v(n)> + ! |uu> is |u(n)> + + ! compute beta(n+1) + beta_i=dznrm2(sze,vv,1) + if (lanczos_debug_print) then + print*,'uu,vv in' + do i=1,n_lanczos_debug + write(6,'(4(E25.15))')uu(i),vv(i) + enddo + endif + ! |vv> is now |u(n+1)> + call zdscal(sze,(1.d0/beta_i),vv,1) + + ! |w(n+1)> = H|u(n+1)> + ! |work> is now |w(n+1)> + call compute_hu(vv,work,sze) + + if (lanczos_debug_print) then + print*,'vv,work' + do i=1,n_lanczos_debug + write(6,'(4(E25.15))')vv(i),work(i) + enddo + endif + + ! alpha(n+1) = + alpha_i=real(u_dot_v_complex(vv,work,sze)) + + do i=1,sze + tmp_c16=work(i)-alpha_i*vv(i)-beta_i*uu(i) + uu(i)=vv(i) + vv(i)=tmp_c16 + enddo + ! |vv> is |v(n+1)> + ! |uu> is |u(n+1)> +end + + + +subroutine lanczos_h(n_lanc_iter,alpha,beta,u1) + implicit none + integer, intent(in) :: n_lanc_iter + double precision, intent(out) :: alpha(n_lanc_iter), beta(n_lanc_iter) + complex*16, intent(in) :: u1(N_det) + integer :: h_size + double precision :: beta_norm, beta_norm_inv + complex*16, allocatable :: vec1(:), vec2(:), vec3(:) + complex*16 :: vec_tmp + double precision, external :: dznrm2 + complex*16, external :: u_dot_v_complex + + integer :: i,j,l + h_size=N_det + BEGIN_DOC + ! lanczos tridiagonalization of H + ! n_lanc_iter is number of lanczos iterations + ! u1 is initial lanczos vector + ! u1 should be normalized + END_DOC + + print *,'starting lanczos' + print *,'h_size = ',h_size +! print *,'initial vector:' +! do i=1,h_size +! print *,u1(i) +! enddo + ! exit if u1 is not normalized + beta_norm = dznrm2(h_size,u1,1) + if (dabs(beta_norm-1.d0) .gt. 1.d-6) then + print *, 'Error: initial Lanczos vector is not normalized' + stop -1 + endif + + allocate(vec1(h_size), & + vec2(h_size), & + vec3(h_size)) + + do i=1,h_size + vec1(i)=u1(i) + enddo + + ! |w1> = H|u1> + ! |vec2> = H|vec1> + call compute_hu(vec1,vec2,h_size)!! TODO: not implemented + + ! alpha(1) = = + ! = + alpha(1)=real(u_dot_v_complex(vec1,vec2,h_size)) + + ! |v1> = |w1> - alpha(1)*|u1> + ! |vec3> = |vec2> - alpha(1)*|vec1> + do i=1,h_size + vec3(i)=vec2(i)-alpha(1)*vec1(i) + enddo + do j=2,n_lanc_iter + call write_time(6) + print *,'starting lanczos iteration:',j + !! vec1 is |u(j-1)> + !! vec3 is |v(j-1)> + + ! beta(j) = sqrt() + beta_norm=dznrm2(h_size,vec3,1) + + ! TODO: check for beta=0? + beta_norm_inv=1.d0/beta_norm + + ! normalize |v(j-1)> to form |u(j)> + call zdscal(h_size,beta_norm_inv,vec3,1) + !! vec3 is |u(j)> + + ! |w(j)> = H|u(j)> + call compute_hu(vec3,vec2,h_size)!! TODO: not implemented + !! vec2 is |w(j)> + + alpha(j)=real(u_dot_v_complex(vec2,vec3,h_size)) + beta(j)=beta_norm + + ! |v(j)> = |w(j)> - alpha(j)*|u(j)> - beta(j)*|u(j-1)> + do l=1,h_size + vec_tmp=vec2(l)-alpha(j)*vec3(l)-beta(j)*vec1(l) + vec1(l)=vec3(l) + vec3(l)=vec_tmp + enddo + !! vec1 is |u(j)> + !! vec3 is |v(j)> + enddo + +end + + +subroutine compute_hu_hp(vec1,vec2,n_hp,h_size,spin_hp,sign_hp,idx_hp) + implicit none + integer, intent(in) :: h_size,n_hp + complex*16, intent(in) :: vec1(h_size,n_hp) + complex*16, intent(out) :: vec2(h_size,n_hp) + integer, intent(in) :: spin_hp(n_hp), idx_hp(n_hp) + double precision, intent (in) :: sign_hp(n_hp) + complex*16 :: vec1_tmp(h_size,n_hp) + integer :: i,j + BEGIN_DOC + ! |vec2> = H|vec1> + ! + ! TODO: implement + ! maybe reuse parts of H_S2_u_0_nstates_{openmp,zmq}? + END_DOC + + vec1_tmp(1:h_size,1:n_hp) = vec1(1:h_size,1:n_hp) + call h_u_0_hp_openmp(vec2,vec1_tmp,n_hp,h_size,spin_hp,sign_hp,idx_hp) + + do j=1,n_hp + do i=1,h_size + if (cdabs(vec1_tmp(i,j) - vec1(i,j)).gt.1.d-6) then + print*,'ERROR: vec1 was changed by h_u_0_openmp' + endif + enddo + enddo +end + +subroutine compute_hu(vec1,vec2,h_size) + implicit none + integer, intent(in) :: h_size + complex*16, intent(in) :: vec1(h_size) + complex*16, intent(out) :: vec2(h_size) + complex*16 :: vec1_tmp(h_size) + integer :: i + BEGIN_DOC + ! |vec2> = H|vec1> + ! + ! TODO: implement + ! maybe reuse parts of H_S2_u_0_nstates_{openmp,zmq}? + END_DOC + + vec1_tmp(1:h_size) = vec1(1:h_size) + call h_u_0_openmp(vec2,vec1_tmp,h_size) + + do i=1,h_size + if (cdabs(vec1_tmp(i) - vec1(i)).gt.1.d-6) then + print*,'ERROR: vec1 was changed by h_u_0_openmp' + endif + enddo +end + +subroutine compute_hu2(vec1,vec2,h_size) + implicit none + integer, intent(in) :: h_size + complex*16, intent(in) :: vec1(h_size) + complex*16, intent(out) :: vec2(h_size) + complex*16, allocatable :: u_tmp(:,:), s_tmp(:,:),v_tmp(:,:) + integer :: i + BEGIN_DOC + ! |vec2> = H|vec1> + ! + ! TODO: implement + ! maybe reuse parts of H_S2_u_0_nstates_{openmp,zmq}? + END_DOC + + allocate(u_tmp(1,h_size),s_tmp(1,h_size),v_tmp(1,h_size)) + + u_tmp(1,1:h_size) = vec1(1:h_size) + call h_s2_u_0_nstates_openmp(v_tmp,s_tmp,u_tmp,1,h_size) + + do i=1,h_size + if (cdabs(u_tmp(1,i) - vec1(i)).gt.1.d-6) then + print*,'ERROR: vec1 was changed by h_u_0_openmp' + endif + enddo + vec2(1:h_size)=v_tmp(1,1:h_size) + deallocate(u_tmp,v_tmp,s_tmp) +end + + + +subroutine diag_lanczos_vals_vecs(alpha, beta, nlanc, vals, vecs, sze) + implicit none + BEGIN_DOC + ! diagonalization of tridiagonal form of H + ! this returns eigenvalues and eigenvectors in vals,vecs + END_DOC + integer, intent(in) :: nlanc,sze + double precision, intent(in) :: alpha(sze), beta(sze) + double precision, intent(out) :: vals(sze), vecs(sze,sze) + double precision :: work(2*nlanc-2), beta_tmp(nlanc-1) + integer :: i,info + + vals(1)=alpha(1) + do i=2,nlanc + vals(i)=alpha(i) + beta_tmp(i-1)=beta(i) + enddo + + call dstev('V', nlanc, vals, beta_tmp, vecs, sze, work, info) + if (info.gt.0) then + print *,'WARNING: diagonalization of tridiagonal form of H did not converge' + else if (info.lt.0) then + print *,'WARNING: argument to dstev had illegal value' + endif +end + +subroutine diag_lanczos_vals_hp(alpha, beta, nlanc, vals, sze,ng) + implicit none + BEGIN_DOC + ! diagonalization of tridiagonal form of H + ! this returns eigenvalues in vals + END_DOC + integer, intent(in) :: nlanc,sze,ng + !double precision, intent(in) :: alpha(ng,sze), beta(sze) + double precision, intent(in) :: alpha(sze,ng), beta(sze,ng) + double precision, intent(out) :: vals(sze,ng) + double precision :: work(1), beta_tmp(nlanc-1), vecs(1) + integer :: i,info,ig + + do ig=1,ng + vals(1,ig)=alpha(1,ig) + do i=2,nlanc + vals(i,ig)=alpha(i,ig) + beta_tmp(i-1)=beta(i,ig) + enddo + + call dstev('N', nlanc, vals(:,ig), beta_tmp, vecs, 1, work, info) + if (info.gt.0) then + print *,'WARNING: diagonalization of tridiagonal form of H did not converge' + else if (info.lt.0) then + print *,'WARNING: argument to dstev had illegal value' + endif + enddo +end +subroutine diag_lanczos_vals(alpha, beta, nlanc, vals, sze) + implicit none + BEGIN_DOC + ! diagonalization of tridiagonal form of H + ! this returns eigenvalues in vals + END_DOC + integer, intent(in) :: nlanc,sze + double precision, intent(in) :: alpha(sze), beta(sze) + double precision, intent(out) :: vals(sze) + double precision :: work(1), beta_tmp(nlanc-1), vecs(1) + integer :: i,info + + vals(1)=alpha(1) + do i=2,nlanc + vals(i)=alpha(i) + beta_tmp(i-1)=beta(i) + enddo + + call dstev('N', nlanc, vals, beta_tmp, vecs, 1, work, info) + if (info.gt.0) then + print *,'WARNING: diagonalization of tridiagonal form of H did not converge' + else if (info.lt.0) then + print *,'WARNING: argument to dstev had illegal value' + endif +end diff --git a/src/green/plot-spec-dens.py b/src/green/plot-spec-dens.py new file mode 100755 index 00000000..88e2dfec --- /dev/null +++ b/src/green/plot-spec-dens.py @@ -0,0 +1,90 @@ +#!/bin/env python + +import gzip +import sys +from math import pi +inv_pi = 1.0/pi + +def spec_dens(alpha,beta,z0,g_sign,e_shift): + sze=len(alpha) + sze_b=len(beta) + if (sze != sze_b): + print('Error: size mismatch',sze,sze_b) + sys.exit(1) + z=z0-g_sign*e_shift + tmp=0.0+0.0j + #for ai,bi in zip(reversed(a),reversed(b)) + for i in range(sze-1,0,-1): + tmp=-(beta[i]**2)/(z+g_sign*alpha[i]+tmp) + tmp=1.0/(z+g_sign*alpha[0]+tmp) + return -1.0 * tmp.imag * inv_pi + +def printspec(ezdir,wmin,wmax,nw,eps): + gdir=ezdir+'/green/' + with open(gdir+'n_green_vec') as infile: + ngvec=int(infile.readline().strip()) + with open(ezdir+'/full_ci_zmq/energy') as infile: + e0=float(infile.readline().strip()) + with open(gdir+'n_lanczos_complete') as infile: + nlanc=int(infile.readline().strip()) + + with gzip.open(gdir+'green_sign.gz') as infile: + gsign0=infile.read().split() + + with gzip.open(gdir+'alpha_lanczos.gz') as infile: + adata0=infile.read().split() + with gzip.open(gdir+'beta_lanczos.gz') as infile: + bdata0=infile.read().split() + + adim=int(adata0.pop(0)) + bdim=int(bdata0.pop(0)) + gsigndim=int(gsign0.pop(0)) + assert adim==2, 'dimension of alpha_lanczos should be 2' + assert bdim==2, 'dimension of beta_lanczos should be 2' + assert gsigndim==1, 'dimension of green_sign should be 1' + + ngvec_2=int(gsign0.pop(0)) + assert ngvec_2==ngvec, 'problem with size of green_sign.gz' + + ashape=tuple(map(int,adata0[:adim])) + bshape=tuple(map(int,bdata0[:bdim])) + assert ashape==(nlanc,ngvec), 'shape of alpha_lanczos should be (nlanc, ngvec)' + assert bshape==(nlanc,ngvec), 'shape of beta_lanczos should be (nlanc, ngvec)' + + amat=[] + for xi in range(ngvec): + amat.append(list(map(float,adata0[adim+xi*nlanc:adim+(xi+1)*nlanc]))) + + bmat=[] + b2mat=[] + for xi in range(ngvec): + #bmat.append(list(map(float,bdata0[bdim+xi*nlanc:bdim+(xi+1)*nlanc]))) + b_tmp=list(map(float,bdata0[bdim+xi*nlanc:bdim+(xi+1)*nlanc])) + b2_tmp=[i*i for i in b_tmp] + bmat.append(b_tmp) + b2mat.append(b2_tmp) + + gsign=list(map(float,gsign0)) + dw=(wmax-wmin)/(nw-1) + wlist = [wmin+iw*dw for iw in range(nw)] + densmat=[] + for ivec in range(ngvec): + densmat.append([spec_dens(amat[ivec],bmat[ivec],iw+1.j*eps,gsign[ivec],e0) for iw in wlist]) + + for i,dd in enumerate(zip(*densmat)): + print(('{:15.6E}'+ngvec*'{:25.15E}').format(wlist[i],*dd)) + +if __name__ == '__main__': + + if len(sys.argv) != 6: + print('bad args') + print('USAGE: plot-spec-dens.py ezfio omega_min omega_max n_omega epsilon') + sys.exit(1) + ezfio=sys.argv[1] + wmin=float(sys.argv[2]) + wmax=float(sys.argv[3]) + nw=int(sys.argv[4]) + eps=float(sys.argv[5]) + printspec(ezfio,wmin,wmax,nw,eps) + + diff --git a/src/green/print_dets_test.irp.f b/src/green/print_dets_test.irp.f new file mode 100644 index 00000000..6466141e --- /dev/null +++ b/src/green/print_dets_test.irp.f @@ -0,0 +1,15 @@ +program print_dets_test + implicit none + read_wf = .True. + touch read_wf + call routine + +end + +subroutine routine + use bitmasks + implicit none + integer :: i + read*,i + print*,psi_det(:,:,i) +end diff --git a/src/green/print_e_mo_debug.irp.f b/src/green/print_e_mo_debug.irp.f new file mode 100644 index 00000000..7bd738bc --- /dev/null +++ b/src/green/print_e_mo_debug.irp.f @@ -0,0 +1,15 @@ +program print_e_mo_debug + implicit none + read_wf = .True. + touch read_wf + call routine + +end + +subroutine routine + use bitmasks + implicit none + integer :: i + read*,i + call print_mo_energies(psi_det(:,:,i),N_int,mo_tot_num) +end diff --git a/src/green/print_h_debug.irp.f b/src/green/print_h_debug.irp.f new file mode 100644 index 00000000..10cc31d3 --- /dev/null +++ b/src/green/print_h_debug.irp.f @@ -0,0 +1,178 @@ +program print_h_debug + implicit none + read_wf = .True. + touch read_wf + call routine + +end + +subroutine routine + use bitmasks + implicit none + integer :: i,j + integer, allocatable :: H_matrix_degree(:,:) + double precision, allocatable :: H_matrix_phase(:,:) + integer :: degree + integer(bit_kind), allocatable :: keys_tmp(:,:,:) + allocate(keys_tmp(N_int,2,N_det)) + do i = 1, N_det + print*,'' + call debug_det(psi_det(1,1,i),N_int) + do j = 1, N_int + keys_tmp(j,1,i) = psi_det(j,1,i) + keys_tmp(j,2,i) = psi_det(j,2,i) + enddo + enddo + if(N_det.gt.10000)then + print*,'Warning !!!' + print*,'Number of determinants is ',N_det + print*,'It means that the H matrix will be enormous !' + print*,'stoppping ..' + stop + endif + print*,'' + print*,'Determinants ' + do i = 1, N_det + enddo + allocate(H_matrix_degree(N_det,N_det),H_matrix_phase(N_det,N_det)) + integer :: exc(0:2,2,2) + double precision :: phase + do i = 1, N_det + do j = i, N_det + call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int) + H_matrix_degree(i,j) = degree + H_matrix_degree(j,i) = degree + phase = 0.d0 + if(degree==1.or.degree==2)then + call get_excitation(psi_det(1,1,i),psi_det(1,1,j),exc,degree,phase,N_int) + endif + H_matrix_phase(i,j) = phase + H_matrix_phase(j,i) = phase + enddo + enddo + print*,'H matrix ' + double precision :: s2 + complex*16 :: ref_h_matrix + ref_h_matrix = h_matrix_all_dets(1,1) + print*,'HF like determinant energy = ',ref_bitmask_energy+nuclear_repulsion + print*,'Ref element of H_matrix = ',ref_h_matrix+nuclear_repulsion + print*,'Printing the H matrix ...' + print*,'' + print*,'' +!do i = 1, N_det +! H_matrix_all_dets(i,i) -= ref_h_matrix +!enddo + + do i = 1, N_det + H_matrix_all_dets(i,i) += nuclear_repulsion + enddo + +!do i = 5,N_det +! H_matrix_all_dets(i,3) = 0.d0 +! H_matrix_all_dets(3,i) = 0.d0 +! H_matrix_all_dets(i,4) = 0.d0 +! H_matrix_all_dets(4,i) = 0.d0 +!enddo + + + + +! TODO: change for complex + do i = 1, N_det + write(*,'(I3,X,A3,2000(E24.15))')i,' | ',H_matrix_all_dets(i,:) + enddo + +! print*,'' +! print*,'' +! print*,'' +! print*,'Printing the degree of excitations within the H matrix' +! print*,'' +! print*,'' +! do i = 1, N_det +! write(*,'(I3,X,A3,X,1000(I1,X))')i,' | ',H_matrix_degree(i,:) +! enddo +! +! +! print*,'' +! print*,'' +! print*,'Printing the phase of the Hamiltonian matrix elements ' +! print*,'' +! print*,'' +! do i = 1, N_det +! write(*,'(I3,X,A3,X,1000(F3.0,X))')i,' | ',H_matrix_phase(i,:) +! enddo +! print*,'' + + +! double precision, allocatable :: eigenvalues(:) +! complex*16, allocatable :: eigenvectors(:,:) +! double precision, allocatable :: s2_eigvalues(:) +! allocate (eigenvectors(size(H_matrix_all_dets,1),N_det)) +! allocate (eigenvalues(N_det),s2_eigvalues(N_det)) +! call lapack_diag_complex(eigenvalues,eigenvectors, & +! H_matrix_all_dets,size(H_matrix_all_dets,1),N_det) +! print*,'Two first eigenvectors ' +! call u_0_S2_u_0(s2_eigvalues,eigenvectors,n_det,keys_tmp,N_int,N_det,size(eigenvectors,1)) +! do j =1, N_states +! print*,'s2 = ',s2_eigvalues(j) +! print*,'e = ',eigenvalues(j) +! print*,'coefs : ' +! do i = 1, N_det +! print*,'i = ',i,eigenvectors(i,j) +! enddo +! if(j>1)then +! print*,'Delta E(H) = ',eigenvalues(1) - eigenvalues(j) +! print*,'Delta E(eV) = ',(eigenvalues(1) - eigenvalues(j))*27.2114d0 +! endif +! enddo +! complex*16 :: get_mo_bielec_integral,k_a_iv,k_b_iv +! integer :: h1,p1,h2,p2 +! h1 = 10 +! p1 = 16 +! h2 = 14 +! p2 = 14 +!!h1 = 1 +!!p1 = 4 +!!h2 = 2 +!!p2 = 2 +! k_a_iv = get_mo_bielec_integral(h1,h2,p2,p1,mo_integrals_map) +! h2 = 15 +! p2 = 15 +! k_b_iv = get_mo_bielec_integral(h1,h2,p2,p1,mo_integrals_map) +! print*,'k_a_iv = ',k_a_iv +! print*,'k_b_iv = ',k_b_iv +! complex*16 :: k_av,k_bv,k_ai,k_bi +! h1 = 16 +! p1 = 14 +! h2 = 14 +! p2 = 16 +! k_av = get_mo_bielec_integral(h1,h2,p1,p2,mo_integrals_map) +! h1 = 16 +! p1 = 15 +! h2 = 15 +! p2 = 16 +! k_bv = get_mo_bielec_integral(h1,h2,p1,p2,mo_integrals_map) +! +! h1 = 10 +! p1 = 14 +! h2 = 14 +! p2 = 10 +! k_ai = get_mo_bielec_integral(h1,h2,p1,p2,mo_integrals_map) +! +! h1 = 10 +! p1 = 15 +! h2 = 15 +! p2 = 10 +! k_bi = get_mo_bielec_integral(h1,h2,p1,p2,mo_integrals_map) +! +! print*,'k_av, k_bv = ',k_av,k_bv +! print*,'k_ai, k_bi = ',k_ai,k_bi +! complex*16 :: k_iv +! +! h1 = 10 +! p1 = 16 +! h2 = 16 +! p2 = 10 +! k_iv = get_mo_bielec_integral(h1,h2,p1,p2,mo_integrals_map) +! print*,'k_iv = ',k_iv +end diff --git a/src/green/print_h_omp_debug.irp.f b/src/green/print_h_omp_debug.irp.f new file mode 100644 index 00000000..abb8b127 --- /dev/null +++ b/src/green/print_h_omp_debug.irp.f @@ -0,0 +1,41 @@ +program print_h_omp_debug + implicit none + read_wf = .True. + touch read_wf + call routine_omp + +end + +subroutine routine_omp + use bitmasks + implicit none + integer :: h_size + complex*16, allocatable :: u_tmp(:,:), s_tmp(:,:),v_tmp(:,:) + integer :: i,n_st + h_size=N_det + BEGIN_DOC + ! |vec2> = H|vec1> + ! + ! TODO: implement + ! maybe reuse parts of H_S2_u_0_nstates_{openmp,zmq}? + END_DOC + n_st=min(1000,h_size) + allocate(u_tmp(n_st,h_size),s_tmp(n_st,h_size),v_tmp(n_st,h_size)) + + u_tmp=(0.d0,0.d0) + v_tmp=(0.d0,0.d0) + s_tmp=(0.d0,0.d0) + + do i=1,n_st + u_tmp(i,i)=(1.d0,0.d0) + enddo + + call h_s2_u_0_nstates_openmp(v_tmp,s_tmp,u_tmp,n_st,h_size) + do i = 1, n_st + v_tmp(i,i) += nuclear_repulsion + enddo + do i = 1, n_st + write(*,'(I3,X,A3,2000(E24.15))')i,' | ',v_tmp(i,:) + enddo + deallocate(u_tmp,v_tmp,s_tmp) +end diff --git a/src/green/print_spectral_dens.irp.f b/src/green/print_spectral_dens.irp.f new file mode 100644 index 00000000..ca6826ce --- /dev/null +++ b/src/green/print_spectral_dens.irp.f @@ -0,0 +1,43 @@ +program print_spectral_dens + implicit none + BEGIN_DOC +! TODO + END_DOC + read_wf = .True. + touch read_wf + provide n_green_vec + call print_lanczos_eigvals + call print_spec +end + +subroutine print_lanczos_eigvals + implicit none + integer :: i, iunit, j + integer :: getunitandopen + character(5) :: jstr + + do j=1,n_green_vec + write(jstr,'(I0.3)') j + iunit = getunitandopen('lanczos_eigval_alpha_beta.out.'//trim(jstr),'w') + print *, 'printing lanczos eigenvalues, alpha, beta to "lanczos_eigval_alpha_beta.out.'//trim(jstr)//'"' + do i=1,n_lanczos_iter + write(iunit,'(I6,3(E25.15))') i, lanczos_eigvals(i,j), alpha_lanczos(i,j), beta_lanczos(i,j) + enddo + close(iunit) + enddo +end +subroutine print_spec + implicit none + integer :: i, iunit, j + integer :: getunitandopen + character(5) :: jstr + do j=1,n_green_vec + write(jstr,'(I0.3)') j + iunit = getunitandopen('omega_A.out.'//trim(jstr),'w') + print *, 'printing frequency, spectral density to "omega_A.out.'//trim(jstr)//'"' + do i=1,n_omega + write(iunit,'(2(E25.15))') omega_list(i), spectral_lanczos(i,j) + enddo + close(iunit) + enddo +end diff --git a/src/green/utils_hp.irp.f b/src/green/utils_hp.irp.f new file mode 100644 index 00000000..0978f9ee --- /dev/null +++ b/src/green/utils_hp.irp.f @@ -0,0 +1,614 @@ +subroutine print_mo_energies(key_ref,nint,nmo) + use bitmasks + BEGIN_DOC + ! get mo energies for one det + END_DOC + implicit none + integer, intent(in) :: nint, nmo + integer(bit_kind), intent(in) :: key_ref(nint,2) + double precision, allocatable :: e_mo(:,:) + integer, allocatable :: occ(:,:),virt(:,:) !(nint*bit_kind_size,2) + integer :: n_occ(2), n_virt(2) + integer, parameter :: int_spin2(1:2) = (/2,1/) + integer :: i,j,ispin,jspin,i0,j0,k + integer(bit_kind), allocatable :: key_virt(:,:) + integer, allocatable :: is_occ(:,:) + + + allocate(occ(nint*bit_kind_size,2),virt(nint*bit_kind_size,2),key_virt(nint,2),e_mo(nmo,2),is_occ(nmo,2)) + is_occ=0 + + call bitstring_to_list_ab(key_ref,occ,n_occ,nint) + do i=1,nint + do ispin=1,2 + key_virt(i,ispin)=xor(full_ijkl_bitmask(i),key_ref(i,ispin)) + enddo + enddo + call bitstring_to_list_ab(key_virt,virt,n_virt,nint) + + e_mo(1:nmo,1)=mo_mono_elec_integral_diag(1:nmo) + e_mo(1:nmo,2)=mo_mono_elec_integral_diag(1:nmo) + + do ispin=1,2 + jspin=int_spin2(ispin) + do i0=1,n_occ(ispin) + i=occ(i0,ispin) + is_occ(i,ispin)=1 + do j0=i0+1,n_occ(ispin) + j=occ(j0,ispin) + e_mo(i,ispin) = e_mo(i,ispin) + mo_bielec_integral_jj_anti(i,j) + e_mo(j,ispin) = e_mo(j,ispin) + mo_bielec_integral_jj_anti(i,j) + enddo + do k=2,ispin + do j0=1,n_occ(jspin) + j=occ(j0,jspin) + e_mo(i,ispin) = e_mo(i,ispin) + mo_bielec_integral_jj(i,j) + e_mo(j,jspin) = e_mo(j,jspin) + mo_bielec_integral_jj(i,j) !can delete this and remove k level of loop + enddo + enddo + do j0=1,n_virt(ispin) + j=virt(j0,ispin) + e_mo(j,ispin) = e_mo(j,ispin) + mo_bielec_integral_jj_anti(i,j) + enddo + do j0=1,n_virt(jspin) + j=virt(j0,jspin) + e_mo(j,jspin) = e_mo(j,jspin) + mo_bielec_integral_jj(i,j) + enddo + enddo + enddo + + do i=1,nmo + write(6,'(2(I5),2(E25.15))')is_occ(i,1),is_occ(i,2),e_mo(i,1),e_mo(i,2) + enddo + deallocate(occ,virt,key_virt,e_mo,is_occ) +end + +subroutine get_mo_energies(key_ref,nint,nmo,e_mo) + use bitmasks + BEGIN_DOC + ! get mo energies for one det + END_DOC + implicit none + integer, intent(in) :: nint, nmo + integer(bit_kind), intent(in) :: key_ref(nint,2) + double precision, intent(out) :: e_mo(nmo,2) + integer, allocatable :: occ(:,:),virt(:,:) !(nint*bit_kind_size,2) + integer :: n_occ(2), n_virt(2) + integer, parameter :: int_spin2(1:2) = (/2,1/) + integer :: i,j,ispin,jspin,i0,j0,k + integer(bit_kind), allocatable :: key_virt(:,:) + + + allocate(occ(nint*bit_kind_size,2),virt(nint*bit_kind_size,2),key_virt(nint,2)) + + call bitstring_to_list_ab(key_ref,occ,n_occ,nint) + do i=1,nint + do ispin=1,2 + key_virt(i,ispin)=xor(full_ijkl_bitmask(i),key_ref(i,ispin)) + enddo + enddo + call bitstring_to_list_ab(key_virt,virt,n_virt,nint) + + e_mo(1:nmo,1)=mo_mono_elec_integral_diag(1:nmo) + e_mo(1:nmo,2)=mo_mono_elec_integral_diag(1:nmo) + + do ispin=1,2 + jspin=int_spin2(ispin) + do i0=1,n_occ(ispin) + i=occ(i0,ispin) + do j0=i0+1,n_occ(ispin) + j=occ(j0,ispin) + e_mo(i,ispin) = e_mo(i,ispin) + mo_bielec_integral_jj_anti(i,j) + e_mo(j,ispin) = e_mo(j,ispin) + mo_bielec_integral_jj_anti(i,j) + enddo + do k=2,ispin + do j0=1,n_occ(jspin) + j=occ(j0,jspin) + e_mo(i,ispin) = e_mo(i,ispin) + mo_bielec_integral_jj(i,j) + e_mo(j,jspin) = e_mo(j,jspin) + mo_bielec_integral_jj(i,j) !can delete this and remove k level of loop + enddo + enddo + do j0=1,n_virt(ispin) + j=virt(j0,ispin) + e_mo(j,ispin) = e_mo(j,ispin) + mo_bielec_integral_jj_anti(i,j) + enddo + do j0=1,n_virt(jspin) + j=virt(j0,jspin) + e_mo(j,jspin) = e_mo(j,jspin) + mo_bielec_integral_jj(i,j) + enddo + enddo + enddo + + deallocate(occ,virt,key_virt) +end + +subroutine get_mask_phase_new(det1, pm, Nint) + use bitmasks + BEGIN_DOC + ! phasemask copied from qp2 + ! return phasemask of det1 in pm + END_DOC + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: det1(Nint,2) + integer(bit_kind), intent(out) :: pm(Nint,2) + integer(bit_kind) :: tmp1, tmp2 + integer :: i + pm(1:Nint,1:2) = det1(1:Nint,1:2) + tmp1 = 0_8 + tmp2 = 0_8 + do i=1,Nint + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 1)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 1)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 2)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 2)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 4)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 4)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 8)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 8)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 16)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 16)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 32)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 32)) + pm(i,1) = ieor(pm(i,1), tmp1) + pm(i,2) = ieor(pm(i,2), tmp2) + if(iand(popcnt(det1(i,1)), 1) == 1) tmp1 = not(tmp1) + if(iand(popcnt(det1(i,2)), 1) == 1) tmp2 = not(tmp2) + end do +end subroutine + +subroutine get_phase_hp(g_idx_int,g_idx_bit,g_spin,g_sign,det_in,g_det_phase,nint,n_g) + use bitmasks + implicit none + integer, intent(in) :: nint,n_g + integer, intent(in) :: g_idx_int(n_g), g_idx_bit(n_g),g_spin(n_g) + double precision, intent(in) :: g_sign(n_g) + integer(bit_kind), intent(in) :: det_in(nint,2) + double precision, intent(out) :: g_det_phase(n_g) + + integer(bit_kind) :: tmp_spindet(nint), pm(nint,2) + double precision, parameter :: phase_dble(0:1) = (/1.d0,-1.d0/) + + integer :: i + logical :: is_allowed(n_g), all_banned, is_filled + + all_banned=.True. + do i=1,n_g + tmp_spindet(1:nint) = det_in(1:nint,g_spin(i)) + call spinorb_is_filled_int_bit(tmp_spindet,g_idx_int(i),g_idx_bit(i),nint,is_filled) + is_allowed(i) = (.not.(((g_sign(i)<0).and.(.not.is_filled)).or.((g_sign(i)>0).and.(is_filled)))) + all_banned=(all_banned.and.(.not.is_allowed(i))) + enddo + + if (all_banned) then + g_det_phase(:)=0.d0 + else + call get_mask_phase_new(det_in,pm,nint) + do i=1,n_g + if (is_allowed(i)) then + g_det_phase(i) = phase_dble(popcnt(iand(ibset(0_bit_kind,g_idx_bit(i)),pm(g_idx_int(i),g_spin(i))))) + else + g_det_phase(i)=0.d0 + endif + enddo + endif +end + +subroutine get_homo_lumo(key_ref,nint,nmo,idx_homo_lumo,spin_homo_lumo) + use bitmasks + implicit none + integer, intent(in) :: nint,nmo + integer(bit_kind), intent(in) :: key_ref(nint,2) + integer, intent(out) :: idx_homo_lumo(2), spin_homo_lumo(2) + + double precision, allocatable :: e_mo(:,:) + integer, allocatable :: occ(:,:),virt(:,:) !(nint*bit_kind_size,2) + integer :: n_occ(2), n_virt(2) + integer :: i,i0,ispin + integer(bit_kind), allocatable :: key_virt(:,:) + double precision :: maxocc(2), minvirt(2) + integer :: imaxocc(2), iminvirt(2) + + allocate(e_mo(nmo,2),key_virt(nint,2),occ(nint*bit_kind_size,2),virt(nint*bit_kind_size,2)) + + call get_mo_energies(key_ref,nint,nmo,e_mo) + + !allocate(occ(nint*bit_kind_size,2),virt(nint*bit_kind_size,2)) + + call bitstring_to_list_ab(key_ref,occ,n_occ,nint) + do i=1,nint + do ispin=1,2 + key_virt(i,ispin)=xor(full_ijkl_bitmask(i),key_ref(i,ispin)) + enddo + enddo + call bitstring_to_list_ab(key_virt,virt,n_virt,nint) + + maxocc=-1.d20 !maybe use -1.d0*huge(0.d0)? + minvirt=1.d20 + imaxocc=-1 + iminvirt=-1 + + do ispin=1,2 + do i0=1,n_occ(ispin) + i=occ(i0,ispin) + if (e_mo(i,ispin).gt.maxocc(ispin)) then + maxocc(ispin)=e_mo(i,ispin) + imaxocc(ispin)=i + endif + enddo + do i0=1,n_virt(ispin) + i=virt(i0,ispin) + if (e_mo(i,ispin).lt.minvirt(ispin)) then + minvirt(ispin)=e_mo(i,ispin) + iminvirt(ispin)=i + endif + enddo + enddo + double precision :: e_mo_thresh + e_mo_thresh = 1.d-8 + !these should both just be 2x2 arrays, but performance here doesn't really matter and this is more readable + !if (maxocc(1).ge.maxocc(2)) then + if ((maxocc(2)-maxocc(1)).le.e_mo_thresh) then + spin_homo_lumo(1)=1 + else + spin_homo_lumo(1)=2 + endif + if ((minvirt(1)-minvirt(2)).le.e_mo_thresh) then + spin_homo_lumo(2)=1 + else + spin_homo_lumo(2)=2 + endif + + idx_homo_lumo(1)=imaxocc(spin_homo_lumo(1)) + idx_homo_lumo(2)=iminvirt(spin_homo_lumo(2)) + + deallocate(e_mo,occ,virt,key_virt) + +end + +subroutine get_list_hp_banned_ab(tmp_det,N_hp,exc_is_banned,spin_hp,sign_hp,idx_hp,nint,all_banned) + use bitmasks + implicit none + BEGIN_DOC + ! input determinant tmp_det and list of single holes/particles + ! for each hole/particle, determine whether it is filled/empty in tmp_det + ! return which are disallowed in exc_is_banned + ! if all are banned, set all_banned to true + END_DOC + integer, intent(in) :: N_hp,nint + integer, intent(in) :: spin_hp(N_hp), idx_hp(N_hp) + double precision, intent(in) :: sign_hp(N_hp) + integer(bit_kind), intent(in) :: tmp_det(nint,2) + logical, intent(out) :: exc_is_banned(N_hp) + logical, intent(out) :: all_banned + + integer :: i + logical :: is_filled + + all_banned = .True. + do i=1,N_hp + call orb_is_filled(tmp_det,idx_hp(i),spin_hp(i),nint,is_filled) + if (sign_hp(i).gt.0) then ! particle creation, banned if already filled + exc_is_banned(i) = is_filled + else ! hole creation, banned if already empty + exc_is_banned(i) = (.not.is_filled) + endif + all_banned = (all_banned.and.exc_is_banned(i)) + enddo +end + +subroutine get_list_hp_banned_single_spin(tmp_spindet,N_hp,exc_is_banned,spin_hp,sign_hp,idx_hp,ispin,nint,all_banned) + use bitmasks + implicit none + BEGIN_DOC + ! input spindeterminant tmp_spindet and list of single holes/particles + ! tmp_spindet is only one spin part of a full det, with spin==ispin + ! for each hole/particle, determine whether it is filled/empty in tmp_det + ! return which are disallowed in exc_is_banned + ! if all are banned, set all_banned to true + END_DOC + integer, intent(in) :: N_hp, ispin, nint + integer, intent(in) :: spin_hp(N_hp), idx_hp(N_hp) + double precision, intent(in) :: sign_hp(N_hp) + integer(bit_kind), intent(in) :: tmp_spindet(nint) + logical, intent(out) :: exc_is_banned(N_hp) + logical, intent(out) :: all_banned + + integer :: i + logical :: is_filled + + all_banned = .True. + do i=1,N_hp + if (spin_hp(i).eq.ispin) then + call orb_is_filled_single_spin(tmp_spindet,idx_hp(i),nint,is_filled) + if (sign_hp(i).gt.0) then ! particle creation, banned if already filled + exc_is_banned(i) = is_filled + else ! hole creation, banned if already empty + exc_is_banned(i) = (.not.is_filled) + endif + else + exc_is_banned(i) = .False. + endif + all_banned = (all_banned.and.exc_is_banned(i)) + enddo +end + +subroutine get_list_hp_banned_spin(tmp_det,N_hp,exc_is_banned,spin_hp,sign_hp,idx_hp,ispin,nint,all_banned) + use bitmasks + implicit none + BEGIN_DOC + ! input determinant tmp_det and list of single holes/particles + ! for each hole/particle, determine whether it is filled/empty in tmp_det + ! return which are disallowed in exc_is_banned + ! if all are banned, set all_banned to true + ! only consider tmp_det(1:N_int, ispin) + END_DOC + integer, intent(in) :: N_hp, ispin, nint + integer, intent(in) :: spin_hp(N_hp), idx_hp(N_hp) + double precision, intent(in) :: sign_hp(N_hp) + integer(bit_kind), intent(in) :: tmp_det(nint,2) + logical, intent(out) :: exc_is_banned(N_hp) + logical, intent(out) :: all_banned + + integer(bit_kind) :: spindet(nint) + + integer :: i + logical :: is_filled + spindet(1:nint) = tmp_det(1:nint,ispin) + + all_banned = .True. + do i=1,N_hp + if (spin_hp(i).eq.ispin) then + call orb_is_filled(tmp_det,idx_hp(i),ispin,nint,is_filled) + if (sign_hp(i).gt.0) then ! particle creation, banned if already filled + exc_is_banned(i) = is_filled + else ! hole creation, banned if already empty + exc_is_banned(i) = (.not.is_filled) + endif + else + exc_is_banned(i) = .False. + endif + all_banned = (all_banned.and.exc_is_banned(i)) + enddo +end + + +subroutine spinorb_is_filled_int_bit(key_ref,iorb_int,iorb_bit,Nint,is_filled) + use bitmasks + implicit none + BEGIN_DOC + ! determine whether iorb is filled in key_ref + ! iorb is specified by int and bit locations within the determinant + END_DOC + integer, intent(in) :: iorb_int, iorb_bit, Nint + integer(bit_kind), intent(in) :: key_ref(Nint) + logical, intent(out) :: is_filled + + ASSERT (iorb_int > 0) + ASSERT (iorb_bit >= 0) + ASSERT (Nint > 0) + is_filled = btest(key_ref(iorb_int),iorb_bit) +end + +subroutine orb_is_filled_int_bit(key_ref,iorb_int,iorb_bit,ispin,Nint,is_filled) + use bitmasks + implicit none + BEGIN_DOC + ! todo: not finished + ! determine whether iorb is filled in key_ref + ! iorb is specified by int and bit locations within the determinant + END_DOC + integer, intent(in) :: iorb_int, iorb_bit, ispin, Nint + integer(bit_kind), intent(in) :: key_ref(Nint,2) + logical, intent(out) :: is_filled + + ASSERT (iorb_int > 0) + ASSERT (iorb_bit >= 0) + ASSERT (ispin > 0) + ASSERT (ispin < 3) + ASSERT (Nint > 0) + is_filled = btest(key_ref(iorb_int,ispin),iorb_bit) +! call spinorb_is_filled_int_bit(key_ref(1,ispin),iorb_int,iorb_bit,Nint,is_filled) +end + +subroutine get_orb_int_bit(iorb,iorb_int,iorb_bit) + BEGIN_DOC + ! get int and bit corresponding to orbital index iorb + END_DOC + use bitmasks + implicit none + integer, intent(in) :: iorb + integer, intent(out) :: iorb_int, iorb_bit + ASSERT (iorb > 0) + iorb_int = ishft(iorb-1,-bit_kind_shift)+1 + ASSERT (iorb_int > 0) + iorb_bit = iorb - ishft(iorb_int-1,bit_kind_shift)-1 + ASSERT (iorb_bit >= 0) +end + +subroutine orb_is_filled_single_spin(key_ref,iorb,Nint,is_filled) + use bitmasks + implicit none + BEGIN_DOC + ! determine whether iorb is filled in key_ref + ! key_ref is single alpha or beta determinant + END_DOC + integer, intent(in) :: iorb, Nint + integer(bit_kind), intent(in) :: key_ref(Nint) + logical, intent(out) :: is_filled + + integer :: k,l + + ASSERT (iorb > 0) + ASSERT (Nint > 0) + + ! k is index of the int where iorb is found + ! l is index of the bit where iorb is found + k = ishft(iorb-1,-bit_kind_shift)+1 + ASSERT (k >0) + l = iorb - ishft(k-1,bit_kind_shift)-1 + ASSERT (l >= 0) + is_filled = btest(key_ref(k),l) +end + +subroutine orb_is_filled(key_ref,iorb,ispin,Nint,is_filled) + use bitmasks + implicit none + BEGIN_DOC + ! determine whether iorb, ispin is filled in key_ref + ! key_ref has alpha and beta parts + END_DOC + integer, intent(in) :: iorb, ispin, Nint + integer(bit_kind), intent(in) :: key_ref(Nint,2) + logical, intent(out) :: is_filled + + integer :: k,l + + ASSERT (iorb > 0) + ASSERT (ispin > 0) + ASSERT (ispin < 3) + ASSERT (Nint > 0) + + ! k is index of the int where iorb is found + ! l is index of the bit where iorb is found + k = ishft(iorb-1,-bit_kind_shift)+1 + ASSERT (k >0) + l = iorb - ishft(k-1,bit_kind_shift)-1 + ASSERT (l >= 0) + is_filled = btest(key_ref(k,ispin),l) +end + +subroutine ac_operator_phase(key_new,key_ref,iorb,ispin,Nint,phase) + use bitmasks + implicit none + BEGIN_DOC + ! apply creation operator to key_ref + ! add electron with spin ispin to orbital with index iorb + ! output resulting det and phase in key_new and phase + END_DOC + integer, intent(in) :: iorb, ispin, Nint + integer(bit_kind), intent(in) :: key_ref(Nint,2) + integer(bit_kind), intent(out) :: key_new(Nint,2) + double precision, intent(out) :: phase + + integer :: k,l,i + + double precision, parameter :: p(0:1) = (/ 1.d0, -1.d0 /) + + ASSERT (iorb > 0) + ASSERT (ispin > 0) + ASSERT (ispin < 3) + ASSERT (Nint > 0) + + key_new=key_ref + + ! alpha det is list of Nint 64-bit ints + ! k is index of the int where iorb is found + ! l is index of the bit where iorb is found + k = ishft(iorb-1,-bit_kind_shift)+1 + ASSERT (k >0) + l = iorb - ishft(k-1,bit_kind_shift)-1 + ASSERT (l >= 0) + key_new(k,ispin) = ibset(key_new(k,ispin),l) + + integer(bit_kind) :: parity_filled + + ! I assume here that the ordering is all alpha spinorbs and then all beta spinorbs + ! If we add an alpha electron, parity is not affected by beta part of determinant + ! (only need number of alpha occupied orbs below iorb) + + ! If we add a beta electron, the parity is affected by alpha part + ! (need total number of occupied alpha orbs (all of which come before beta) + ! and total number of beta occupied orbs below iorb) + + if (ispin==1) then + parity_filled=0_bit_kind + else + parity_filled=iand(elec_alpha_num,1_bit_kind) + endif + + ! get parity due to orbs in other ints (with lower indices) + do i=1,k-1 + parity_filled = iand(popcnt(key_ref(i,ispin)),parity_filled) + enddo + + ! get parity due to orbs in same int as iorb + ! ishft(1_bit_kind,l)-1 has its l rightmost bits set to 1, other bits set to 0 + parity_filled = iand(popcnt(iand(ishft(1_bit_kind,l)-1,key_ref(k,ispin))),parity_filled) + phase = p(iand(1_bit_kind,parity_filled)) + +end + +subroutine a_operator_phase(key_new,key_ref,iorb,ispin,Nint,phase) + use bitmasks + implicit none + BEGIN_DOC + ! apply annihilation operator to key_ref + ! remove electron with spin ispin to orbital with index iorb + ! output resulting det and phase in key_new and phase + END_DOC + integer, intent(in) :: iorb, ispin, Nint + integer(bit_kind), intent(in) :: key_ref(Nint,2) + integer(bit_kind), intent(out) :: key_new(Nint,2) + double precision, intent(out) :: phase + + integer :: k,l,i + + double precision, parameter :: p(0:1) = (/ 1.d0, -1.d0 /) + + ASSERT (iorb > 0) + ASSERT (ispin > 0) + ASSERT (ispin < 3) + ASSERT (Nint > 0) + + key_new=key_ref + + ! alpha det is list of Nint 64-bit ints + ! k is index of the int where iorb is found + ! l is index of the bit where iorb is found + k = ishft(iorb-1,-bit_kind_shift)+1 + ASSERT (k >0) + l = iorb - ishft(k-1,bit_kind_shift)-1 + ASSERT (l >= 0) + key_new(k,ispin) = ibclr(key_new(k,ispin),l) + + integer(bit_kind) :: parity_filled + + ! I assume here that the ordering is all alpha spinorbs and then all beta spinorbs + ! If we add an alpha electron, parity is not affected by beta part of determinant + ! (only need number of alpha occupied orbs below iorb) + + ! If we add a beta electron, the parity is affected by alpha part + ! (need total number of occupied alpha orbs (all of which come before beta) + ! and total number of beta occupied orbs below iorb) + + if (ispin==1) then + parity_filled=0_bit_kind + else + parity_filled=iand(elec_alpha_num,1_bit_kind) + endif + + ! get parity due to orbs in other ints (with lower indices) + do i=1,k-1 + parity_filled = iand(popcnt(key_ref(i,ispin)),parity_filled) + enddo + + ! get parity due to orbs in same int as iorb + ! ishft(1_bit_kind,l)-1 has its l rightmost bits set to 1, other bits set to 0 + parity_filled = iand(popcnt(iand(ishft(1_bit_kind,l)-1,key_ref(k,ispin))),parity_filled) + phase = p(iand(1_bit_kind,parity_filled)) + +end +BEGIN_PROVIDER [ double precision, mo_mono_elec_integral_diag,(mo_tot_num)] + implicit none + integer :: i + BEGIN_DOC + ! diagonal elements of mo_mono_elec_integral array + END_DOC + print*,'Providing the mono electronic integrals (diagonal)' + + do i = 1, mo_tot_num + mo_mono_elec_integral_diag(i) = real(mo_mono_elec_integral(i,i)) + enddo + +END_PROVIDER From 25d0cbaa7537647c08a02fd51a61ec6c37116dca Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Tue, 2 Jun 2020 11:59:14 -0500 Subject: [PATCH 2/8] complex to double in ezfio --- src/green/EZFIO.cfg | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/green/EZFIO.cfg b/src/green/EZFIO.cfg index 859a3eb9..0ad57888 100644 --- a/src/green/EZFIO.cfg +++ b/src/green/EZFIO.cfg @@ -73,14 +73,14 @@ size: (green.n_lanczos_iter,green.n_green_vec) [un_lanczos] interface: ezfio doc: saved lanczos u vector -type: complex*16 -size: (determinants.n_det,green.n_green_vec) +type: double precision +size: (2,determinants.n_det,green.n_green_vec) [vn_lanczos] interface: ezfio doc: saved lanczos v vector -type: complex*16 -size: (determinants.n_det,green.n_green_vec) +type: double precision +size: (2,determinants.n_det,green.n_green_vec) [lanczos_eigvals] interface: ezfio From 0fd6eb3897a9abff08574b0f58eb3d4ca6d35784 Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Wed, 3 Jun 2020 16:13:16 -0500 Subject: [PATCH 3/8] updated green for qp2 --- src/green/green.main.irp.f | 2 +- src/green/hu0_hp.irp.f | 54 ++++++++++++------------ src/green/hu0_lanczos.irp.f | 10 ++--- src/green/lanczos.irp.f | 13 +++--- src/green/plot-spec-dens.py | 2 +- src/green/print_e_mo_debug.irp.f | 2 +- src/green/print_h_debug.irp.f | 6 +-- src/green/print_h_omp_debug.irp.f | 2 +- src/green/utils_hp.irp.f | 70 +++++++++++++++---------------- src/utils/constants.include.F | 1 + 10 files changed, 82 insertions(+), 80 deletions(-) diff --git a/src/green/green.main.irp.f b/src/green/green.main.irp.f index c9b3ef66..3fe26424 100644 --- a/src/green/green.main.irp.f +++ b/src/green/green.main.irp.f @@ -15,7 +15,7 @@ end subroutine psicoefprinttest implicit none integer :: i - TOUCH psi_coef + TOUCH psi_coef_complex print *, 'printing ndet', N_det end subroutine print_lanczos_eigvals diff --git a/src/green/hu0_hp.irp.f b/src/green/hu0_hp.irp.f index c3d8be40..4fa7275f 100644 --- a/src/green/hu0_hp.irp.f +++ b/src/green/hu0_hp.irp.f @@ -595,22 +595,22 @@ subroutine i_h_j_double_spin_hp(key_i,key_j,Nint,ispin,hij_hp,N_hp,spin_hp,sign_ double precision :: phase_hp(N_hp) integer :: exc(0:2,2) double precision :: phase - complex*16, external :: get_mo_bielec_integral + complex*16, external :: mo_two_e_integral_complex integer :: i1,i2,i3,i4,j2,j3,ii - PROVIDE big_array_exchange_integrals mo_bielec_integrals_in_map + PROVIDE big_array_exchange_integrals_complex mo_two_e_integrals_in_map call get_double_excitation_spin(key_i,key_j,exc,phase,Nint) - hij0 = phase*(get_mo_bielec_integral( & + hij0 = phase*(mo_two_e_integral_complex( & exc(1,1), & exc(2,1), & exc(1,2), & - exc(2,2), mo_integrals_map) - & - get_mo_bielec_integral( & + exc(2,2)) - & + mo_two_e_integral_complex( & exc(1,1), & exc(2,1), & exc(2,2), & - exc(1,2), mo_integrals_map) ) + exc(1,2)) ) ASSERT (exc(1,1) < exc(2,1)) ASSERT (exc(1,2) < exc(2,2)) @@ -661,14 +661,14 @@ subroutine i_h_j_mono_spin_hp(key_i,key_j,Nint,spin,hij_hp,N_hp,spin_hp,sign_hp, integer :: exc(0:2,2) double precision :: phase - PROVIDE big_array_exchange_integrals mo_bielec_integrals_in_map + PROVIDE big_array_exchange_integrals_complex mo_two_e_integrals_in_map - call get_mono_excitation_spin(key_i(1,spin),key_j(1,spin),exc,phase,Nint) + call get_single_excitation_spin(key_i(1,spin),key_j(1,spin),exc,phase,Nint) - call get_mono_excitation_from_fock_hp(key_i,key_j,exc(1,1),exc(1,2),spin,phase,N_hp,hij_hp,spin_hp,sign_hp,idx_hp,allowed_hp) + call get_single_excitation_from_fock_hp(key_i,key_j,exc(1,1),exc(1,2),spin,phase,N_hp,hij_hp,spin_hp,sign_hp,idx_hp,allowed_hp) end -subroutine get_mono_excitation_from_fock_hp(det_1,det_2,h,p,spin,phase,N_hp,hij_hp,spin_hp,sign_hp,idx_hp,allowed_hp) +subroutine get_single_excitation_from_fock_hp(det_1,det_2,h,p,spin,phase,N_hp,hij_hp,spin_hp,sign_hp,idx_hp,allowed_hp) use bitmasks implicit none integer,intent(in) :: h,p,spin,N_hp @@ -699,37 +699,37 @@ subroutine get_mono_excitation_from_fock_hp(det_1,det_2,h,p,spin,phase,N_hp,hij_ enddo call bitstring_to_list_ab(hole, occ_hole, n_occ_ab_hole, N_int) call bitstring_to_list_ab(partcl, occ_partcl, n_occ_ab_partcl, N_int) - hij0 = fock_operator_closed_shell_ref_bitmask(h,p) + hij0 = fock_op_cshell_ref_bitmask_cplx(h,p) ! holes :: direct terms do i0 = 1, n_occ_ab_hole(1) i = occ_hole(i0,1) - hij0 -= big_array_coulomb_integrals(i,h,p) ! get_mo_bielec_integral_schwartz(h,i,p,i,mo_integrals_map) + hij0 -= big_array_coulomb_integrals_complex(i,h,p) ! get_mo_bielec_integral_schwartz(h,i,p,i,mo_integrals_map) enddo do i0 = 1, n_occ_ab_hole(2) i = occ_hole(i0,2) - hij0 -= big_array_coulomb_integrals(i,h,p) !get_mo_bielec_integral_schwartz(h,i,p,i,mo_integrals_map) + hij0 -= big_array_coulomb_integrals_complex(i,h,p) !get_mo_bielec_integral_schwartz(h,i,p,i,mo_integrals_map) enddo ! holes :: exchange terms do i0 = 1, n_occ_ab_hole(spin) i = occ_hole(i0,spin) - hij0 += big_array_exchange_integrals(i,h,p) ! get_mo_bielec_integral_schwartz(h,i,i,p,mo_integrals_map) + hij0 += big_array_exchange_integrals_complex(i,h,p) ! get_mo_bielec_integral_schwartz(h,i,i,p,mo_integrals_map) enddo ! particles :: direct terms do i0 = 1, n_occ_ab_partcl(1) i = occ_partcl(i0,1) - hij0 += big_array_coulomb_integrals(i,h,p)!get_mo_bielec_integral_schwartz(h,i,p,i,mo_integrals_map) + hij0 += big_array_coulomb_integrals_complex(i,h,p)!get_mo_bielec_integral_schwartz(h,i,p,i,mo_integrals_map) enddo do i0 = 1, n_occ_ab_partcl(2) i = occ_partcl(i0,2) - hij0 += big_array_coulomb_integrals(i,h,p) !get_mo_bielec_integral_schwartz(h,i,p,i,mo_integrals_map) + hij0 += big_array_coulomb_integrals_complex(i,h,p) !get_mo_bielec_integral_schwartz(h,i,p,i,mo_integrals_map) enddo ! particles :: exchange terms do i0 = 1, n_occ_ab_partcl(spin) i = occ_partcl(i0,spin) - hij0 -= big_array_exchange_integrals(i,h,p)!get_mo_bielec_integral_schwartz(h,i,i,p,mo_integrals_map) + hij0 -= big_array_exchange_integrals_complex(i,h,p)!get_mo_bielec_integral_schwartz(h,i,i,p,mo_integrals_map) enddo low=min(h,p) @@ -771,7 +771,7 @@ subroutine get_mono_excitation_from_fock_hp(det_1,det_2,h,p,spin,phase,N_hp,hij_ hij_hp(ii) = 0.d0 cycle else if (spin.eq.spin_hp(ii)) then - hij_hp(ii) = hij0 + sign_hp(ii) *(big_array_coulomb_integrals(idx_hp(ii),h,p) - big_array_exchange_integrals(idx_hp(ii),h,p)) + hij_hp(ii) = hij0 + sign_hp(ii) *(big_array_coulomb_integrals_complex(idx_hp(ii),h,p) - big_array_exchange_integrals_complex(idx_hp(ii),h,p)) if ((low.lt.idx_hp(ii)).and.(high.gt.idx_hp(ii))) then phase_hp(ii) = -1.d0 else @@ -779,7 +779,7 @@ subroutine get_mono_excitation_from_fock_hp(det_1,det_2,h,p,spin,phase,N_hp,hij_ endif else phase_hp(ii) = 1.d0 - hij_hp(ii) = hij0 + sign_hp(ii) * big_array_coulomb_integrals(idx_hp(ii),h,p) + hij_hp(ii) = hij0 + sign_hp(ii) * big_array_coulomb_integrals_complex(idx_hp(ii),h,p) endif hij_hp(ii) = hij_hp(ii) * phase * phase_hp(ii) enddo @@ -806,24 +806,24 @@ subroutine i_H_j_double_alpha_beta_hp(key_i,key_j,Nint,hij_hp,N_hp,spin_hp,sign_ integer :: lowhigh(2,2) integer :: exc(0:2,2,2) double precision :: phase, phase2 - complex*16, external :: get_mo_bielec_integral + complex*16, external :: mo_two_e_integral_complex - PROVIDE big_array_exchange_integrals mo_bielec_integrals_in_map + PROVIDE big_array_exchange_integrals_complex mo_two_e_integrals_in_map - call get_mono_excitation_spin(key_i(1,1),key_j(1,1),exc(0,1,1),phase,Nint) - call get_mono_excitation_spin(key_i(1,2),key_j(1,2),exc(0,1,2),phase2,Nint) + call get_single_excitation_spin(key_i(1,1),key_j(1,1),exc(0,1,1),phase,Nint) + call get_single_excitation_spin(key_i(1,2),key_j(1,2),exc(0,1,2),phase2,Nint) phase = phase*phase2 if (exc(1,1,1) == exc(1,2,2)) then - hij0 = big_array_exchange_integrals(exc(1,1,1),exc(1,1,2),exc(1,2,1)) + hij0 = big_array_exchange_integrals_complex(exc(1,1,1),exc(1,1,2),exc(1,2,1)) else if (exc(1,2,1) == exc(1,1,2)) then - hij0 = big_array_exchange_integrals(exc(1,2,1),exc(1,1,1),exc(1,2,2)) + hij0 = big_array_exchange_integrals_complex(exc(1,2,1),exc(1,1,1),exc(1,2,2)) else - hij0 = get_mo_bielec_integral( & + hij0 = mo_two_e_integral_complex( & exc(1,1,1), & exc(1,1,2), & exc(1,2,1), & - exc(1,2,2) ,mo_integrals_map) + exc(1,2,2)) endif !todo: clean this up diff --git a/src/green/hu0_lanczos.irp.f b/src/green/hu0_lanczos.irp.f index e4da5c78..6f7ebf1d 100644 --- a/src/green/hu0_lanczos.irp.f +++ b/src/green/hu0_lanczos.irp.f @@ -194,7 +194,7 @@ subroutine H_u_0_openmp_work_$N_int(v_t,u_t,sze,istart,iend,ishift,istep) ASSERT (lrow <= N_det_alpha_unique) tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow) - call i_h_j_double_alpha_beta(tmp_det,tmp_det2,$N_int,hij) + call i_h_j_double_alpha_beta_complex(tmp_det,tmp_det2,$N_int,hij) v_t(k_a) = v_t(k_a) + hij * u_t(l_a) enddo enddo @@ -264,7 +264,7 @@ subroutine H_u_0_openmp_work_$N_int(v_t,u_t,sze,istart,iend,ishift,istep) ASSERT (lrow <= N_det_alpha_unique) tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow) - call i_H_j_mono_spin( tmp_det, tmp_det2, $N_int, 1, hij) + call i_h_j_single_spin_complex( tmp_det, tmp_det2, $N_int, 1, hij) v_t(k_a) = v_t(k_a) + hij * u_t(l_a) enddo @@ -280,7 +280,7 @@ subroutine H_u_0_openmp_work_$N_int(v_t,u_t,sze,istart,iend,ishift,istep) lrow = psi_bilinear_matrix_rows(l_a) ASSERT (lrow <= N_det_alpha_unique) - call i_H_j_double_spin( tmp_det(1,1), psi_det_alpha_unique(1, lrow), $N_int, hij) + call i_h_j_double_spin_complex( tmp_det(1,1), psi_det_alpha_unique(1, lrow), $N_int, hij) v_t(k_a) = v_t(k_a) + hij * u_t(l_a) enddo @@ -341,7 +341,7 @@ subroutine H_u_0_openmp_work_$N_int(v_t,u_t,sze,istart,iend,ishift,istep) ASSERT (lcol <= N_det_beta_unique) tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, lcol) - call i_H_j_mono_spin( tmp_det, tmp_det2, $N_int, 2, hij) + call i_h_j_single_spin_complex( tmp_det, tmp_det2, $N_int, 2, hij) l_a = psi_bilinear_matrix_transp_order(l_b) ASSERT (l_a <= N_det) v_t(k_a) = v_t(k_a) + hij * u_t(l_a) @@ -357,7 +357,7 @@ subroutine H_u_0_openmp_work_$N_int(v_t,u_t,sze,istart,iend,ishift,istep) lcol = psi_bilinear_matrix_transp_columns(l_b) ASSERT (lcol <= N_det_beta_unique) - call i_H_j_double_spin( tmp_det(1,2), psi_det_beta_unique(1, lcol), $N_int, hij) + call i_h_j_double_spin_complex( tmp_det(1,2), psi_det_beta_unique(1, lcol), $N_int, hij) l_a = psi_bilinear_matrix_transp_order(l_b) ASSERT (l_a <= N_det) diff --git a/src/green/lanczos.irp.f b/src/green/lanczos.irp.f index a2557abb..baf66d80 100644 --- a/src/green/lanczos.irp.f +++ b/src/green/lanczos.irp.f @@ -28,7 +28,7 @@ END_PROVIDER integer :: idx_homo_lumo(2), spin_homo_lumo(2) logical :: has_idx,has_spin,has_sign,has_lanc integer :: nlanc - ! needs psi_det, mo_tot_num, N_int, mo_bielec_integral_jj, mo_mono_elec_integral_diag + ! needs psi_det, mo_num, N_int, mo_bielec_integral_jj, mo_mono_elec_integral_diag call ezfio_has_green_green_idx(has_idx) call ezfio_has_green_green_spin(has_spin) call ezfio_has_green_green_sign(has_sign) @@ -43,7 +43,7 @@ END_PROVIDER stop 'problem with lanczos restart; need idx, spin, sign' else print*,'new lanczos calculation, finding homo/lumo' - call get_homo_lumo(psi_det(1:N_int,1:2,1),N_int,mo_tot_num,idx_homo_lumo,spin_homo_lumo) + call get_homo_lumo(psi_det(1:N_int,1:2,1),N_int,mo_num,idx_homo_lumo,spin_homo_lumo) ! homo green_idx(1)=idx_homo_lumo(1) @@ -75,7 +75,7 @@ END_PROVIDER ! endif ! else ! print*,'new lanczos calculation, finding homo/lumo' -! call get_homo_lumo(psi_det(1:N_int,1:2,1),N_int,mo_tot_num,idx_homo_lumo,spin_homo_lumo) +! call get_homo_lumo(psi_det(1:N_int,1:2,1),N_int,mo_num,idx_homo_lumo,spin_homo_lumo) ! ! ! homo ! green_idx(1)=idx_homo_lumo(1) @@ -279,9 +279,9 @@ BEGIN_PROVIDER [ double precision, spectral_lanczos, (n_omega,n_green_vec) ] logical :: has_ci_energy double precision :: ref_energy_0 PROVIDE delta_omega alpha_lanczos beta_lanczos omega_list - call ezfio_has_full_ci_zmq_energy(has_ci_energy) + call ezfio_has_fci_energy(has_ci_energy) if (has_ci_energy) then - call ezfio_get_full_ci_zmq_energy(ref_energy_0) + call ezfio_get_fci_energy(ref_energy_0) else print*,'no reference energy from full_ci_zmq, exiting' stop @@ -469,7 +469,8 @@ subroutine lanczos_h_step_hp(uu,vv,work,sze,alpha_i,beta_i,ng,spin_hp,sign_hp,id complex*16, intent(inout) :: uu(sze,ng),vv(sze,ng) complex*16, intent(out) :: work(sze,ng) double precision, intent(out) :: alpha_i(ng), beta_i(ng) - integer, intent(in) :: spin_hp(ng), sign_hp(ng), idx_hp(ng) + integer, intent(in) :: spin_hp(ng), idx_hp(ng) + double precision, intent(in) :: sign_hp(ng) double precision, external :: dznrm2 complex*16, external :: u_dot_v_complex diff --git a/src/green/plot-spec-dens.py b/src/green/plot-spec-dens.py index 88e2dfec..bf4f2294 100755 --- a/src/green/plot-spec-dens.py +++ b/src/green/plot-spec-dens.py @@ -23,7 +23,7 @@ def printspec(ezdir,wmin,wmax,nw,eps): gdir=ezdir+'/green/' with open(gdir+'n_green_vec') as infile: ngvec=int(infile.readline().strip()) - with open(ezdir+'/full_ci_zmq/energy') as infile: + with open(ezdir+'/fci/energy') as infile: e0=float(infile.readline().strip()) with open(gdir+'n_lanczos_complete') as infile: nlanc=int(infile.readline().strip()) diff --git a/src/green/print_e_mo_debug.irp.f b/src/green/print_e_mo_debug.irp.f index 7bd738bc..1fe41e34 100644 --- a/src/green/print_e_mo_debug.irp.f +++ b/src/green/print_e_mo_debug.irp.f @@ -11,5 +11,5 @@ subroutine routine implicit none integer :: i read*,i - call print_mo_energies(psi_det(:,:,i),N_int,mo_tot_num) + call print_mo_energies(psi_det(:,:,i),N_int,mo_num) end diff --git a/src/green/print_h_debug.irp.f b/src/green/print_h_debug.irp.f index 10cc31d3..4dd394d7 100644 --- a/src/green/print_h_debug.irp.f +++ b/src/green/print_h_debug.irp.f @@ -53,7 +53,7 @@ subroutine routine print*,'H matrix ' double precision :: s2 complex*16 :: ref_h_matrix - ref_h_matrix = h_matrix_all_dets(1,1) + ref_h_matrix = h_matrix_all_dets_complex(1,1) print*,'HF like determinant energy = ',ref_bitmask_energy+nuclear_repulsion print*,'Ref element of H_matrix = ',ref_h_matrix+nuclear_repulsion print*,'Printing the H matrix ...' @@ -64,7 +64,7 @@ subroutine routine !enddo do i = 1, N_det - H_matrix_all_dets(i,i) += nuclear_repulsion + H_matrix_all_dets_complex(i,i) += nuclear_repulsion enddo !do i = 5,N_det @@ -79,7 +79,7 @@ subroutine routine ! TODO: change for complex do i = 1, N_det - write(*,'(I3,X,A3,2000(E24.15))')i,' | ',H_matrix_all_dets(i,:) + write(*,'(I3,X,A3,2000(E24.15))')i,' | ',H_matrix_all_dets_complex(i,:) enddo ! print*,'' diff --git a/src/green/print_h_omp_debug.irp.f b/src/green/print_h_omp_debug.irp.f index abb8b127..0a9cd930 100644 --- a/src/green/print_h_omp_debug.irp.f +++ b/src/green/print_h_omp_debug.irp.f @@ -30,7 +30,7 @@ subroutine routine_omp u_tmp(i,i)=(1.d0,0.d0) enddo - call h_s2_u_0_nstates_openmp(v_tmp,s_tmp,u_tmp,n_st,h_size) + call h_s2_u_0_nstates_openmp_complex(v_tmp,s_tmp,u_tmp,n_st,h_size) do i = 1, n_st v_tmp(i,i) += nuclear_repulsion enddo diff --git a/src/green/utils_hp.irp.f b/src/green/utils_hp.irp.f index 0978f9ee..264e3014 100644 --- a/src/green/utils_hp.irp.f +++ b/src/green/utils_hp.irp.f @@ -26,8 +26,8 @@ subroutine print_mo_energies(key_ref,nint,nmo) enddo call bitstring_to_list_ab(key_virt,virt,n_virt,nint) - e_mo(1:nmo,1)=mo_mono_elec_integral_diag(1:nmo) - e_mo(1:nmo,2)=mo_mono_elec_integral_diag(1:nmo) + e_mo(1:nmo,1)=mo_one_e_integrals_diag(1:nmo) + e_mo(1:nmo,2)=mo_one_e_integrals_diag(1:nmo) do ispin=1,2 jspin=int_spin2(ispin) @@ -36,23 +36,23 @@ subroutine print_mo_energies(key_ref,nint,nmo) is_occ(i,ispin)=1 do j0=i0+1,n_occ(ispin) j=occ(j0,ispin) - e_mo(i,ispin) = e_mo(i,ispin) + mo_bielec_integral_jj_anti(i,j) - e_mo(j,ispin) = e_mo(j,ispin) + mo_bielec_integral_jj_anti(i,j) + e_mo(i,ispin) = e_mo(i,ispin) + mo_two_e_integrals_jj_anti(i,j) + e_mo(j,ispin) = e_mo(j,ispin) + mo_two_e_integrals_jj_anti(i,j) enddo do k=2,ispin do j0=1,n_occ(jspin) j=occ(j0,jspin) - e_mo(i,ispin) = e_mo(i,ispin) + mo_bielec_integral_jj(i,j) - e_mo(j,jspin) = e_mo(j,jspin) + mo_bielec_integral_jj(i,j) !can delete this and remove k level of loop + e_mo(i,ispin) = e_mo(i,ispin) + mo_two_e_integrals_jj(i,j) + e_mo(j,jspin) = e_mo(j,jspin) + mo_two_e_integrals_jj(i,j) !can delete this and remove k level of loop enddo enddo do j0=1,n_virt(ispin) j=virt(j0,ispin) - e_mo(j,ispin) = e_mo(j,ispin) + mo_bielec_integral_jj_anti(i,j) + e_mo(j,ispin) = e_mo(j,ispin) + mo_two_e_integrals_jj_anti(i,j) enddo do j0=1,n_virt(jspin) j=virt(j0,jspin) - e_mo(j,jspin) = e_mo(j,jspin) + mo_bielec_integral_jj(i,j) + e_mo(j,jspin) = e_mo(j,jspin) + mo_two_e_integrals_jj(i,j) enddo enddo enddo @@ -89,8 +89,8 @@ subroutine get_mo_energies(key_ref,nint,nmo,e_mo) enddo call bitstring_to_list_ab(key_virt,virt,n_virt,nint) - e_mo(1:nmo,1)=mo_mono_elec_integral_diag(1:nmo) - e_mo(1:nmo,2)=mo_mono_elec_integral_diag(1:nmo) + e_mo(1:nmo,1)=mo_one_e_integrals_diag(1:nmo) + e_mo(1:nmo,2)=mo_one_e_integrals_diag(1:nmo) do ispin=1,2 jspin=int_spin2(ispin) @@ -98,23 +98,23 @@ subroutine get_mo_energies(key_ref,nint,nmo,e_mo) i=occ(i0,ispin) do j0=i0+1,n_occ(ispin) j=occ(j0,ispin) - e_mo(i,ispin) = e_mo(i,ispin) + mo_bielec_integral_jj_anti(i,j) - e_mo(j,ispin) = e_mo(j,ispin) + mo_bielec_integral_jj_anti(i,j) + e_mo(i,ispin) = e_mo(i,ispin) + mo_two_e_integrals_jj_anti(i,j) + e_mo(j,ispin) = e_mo(j,ispin) + mo_two_e_integrals_jj_anti(i,j) enddo do k=2,ispin do j0=1,n_occ(jspin) j=occ(j0,jspin) - e_mo(i,ispin) = e_mo(i,ispin) + mo_bielec_integral_jj(i,j) - e_mo(j,jspin) = e_mo(j,jspin) + mo_bielec_integral_jj(i,j) !can delete this and remove k level of loop + e_mo(i,ispin) = e_mo(i,ispin) + mo_two_e_integrals_jj(i,j) + e_mo(j,jspin) = e_mo(j,jspin) + mo_two_e_integrals_jj(i,j) !can delete this and remove k level of loop enddo enddo do j0=1,n_virt(ispin) j=virt(j0,ispin) - e_mo(j,ispin) = e_mo(j,ispin) + mo_bielec_integral_jj_anti(i,j) + e_mo(j,ispin) = e_mo(j,ispin) + mo_two_e_integrals_jj_anti(i,j) enddo do j0=1,n_virt(jspin) j=virt(j0,jspin) - e_mo(j,jspin) = e_mo(j,jspin) + mo_bielec_integral_jj(i,j) + e_mo(j,jspin) = e_mo(j,jspin) + mo_two_e_integrals_jj(i,j) enddo enddo enddo @@ -524,17 +524,17 @@ subroutine ac_operator_phase(key_new,key_ref,iorb,ispin,Nint,phase) if (ispin==1) then parity_filled=0_bit_kind else - parity_filled=iand(elec_alpha_num,1_bit_kind) + parity_filled=iand(int(elec_alpha_num,bit_kind),1_bit_kind) endif ! get parity due to orbs in other ints (with lower indices) do i=1,k-1 - parity_filled = iand(popcnt(key_ref(i,ispin)),parity_filled) + parity_filled = iand(int(popcnt(key_ref(i,ispin)),bit_kind),parity_filled) enddo ! get parity due to orbs in same int as iorb ! ishft(1_bit_kind,l)-1 has its l rightmost bits set to 1, other bits set to 0 - parity_filled = iand(popcnt(iand(ishft(1_bit_kind,l)-1,key_ref(k,ispin))),parity_filled) + parity_filled = iand(int(popcnt(iand(ishft(1_bit_kind,l)-1,key_ref(k,ispin))),bit_kind),parity_filled) phase = p(iand(1_bit_kind,parity_filled)) end @@ -585,30 +585,30 @@ subroutine a_operator_phase(key_new,key_ref,iorb,ispin,Nint,phase) if (ispin==1) then parity_filled=0_bit_kind else - parity_filled=iand(elec_alpha_num,1_bit_kind) + parity_filled=iand(int(elec_alpha_num,bit_kind),1_bit_kind) endif ! get parity due to orbs in other ints (with lower indices) do i=1,k-1 - parity_filled = iand(popcnt(key_ref(i,ispin)),parity_filled) + parity_filled = iand(int(popcnt(key_ref(i,ispin)),bit_kind),parity_filled) enddo ! get parity due to orbs in same int as iorb ! ishft(1_bit_kind,l)-1 has its l rightmost bits set to 1, other bits set to 0 - parity_filled = iand(popcnt(iand(ishft(1_bit_kind,l)-1,key_ref(k,ispin))),parity_filled) + parity_filled = iand(int(popcnt(iand(ishft(1_bit_kind,l)-1,key_ref(k,ispin))),bit_kind),parity_filled) phase = p(iand(1_bit_kind,parity_filled)) end -BEGIN_PROVIDER [ double precision, mo_mono_elec_integral_diag,(mo_tot_num)] - implicit none - integer :: i - BEGIN_DOC - ! diagonal elements of mo_mono_elec_integral array - END_DOC - print*,'Providing the mono electronic integrals (diagonal)' - - do i = 1, mo_tot_num - mo_mono_elec_integral_diag(i) = real(mo_mono_elec_integral(i,i)) - enddo - -END_PROVIDER +!BEGIN_PROVIDER [ double precision, mo_mono_elec_integral_diag,(mo_num)] +! implicit none +! integer :: i +! BEGIN_DOC +! ! diagonal elements of mo_mono_elec_integral array +! END_DOC +! print*,'Providing the mono electronic integrals (diagonal)' +! +! do i = 1, mo_num +! mo_mono_elec_integral_diag(i) = real(mo_mono_elec_integral(i,i)) +! enddo +! +!END_PROVIDER diff --git a/src/utils/constants.include.F b/src/utils/constants.include.F index 7399b4a6..bad68054 100644 --- a/src/utils/constants.include.F +++ b/src/utils/constants.include.F @@ -7,6 +7,7 @@ double precision, parameter :: sqpi = dsqrt(dacos(-1.d0)) double precision, parameter :: pi_5_2 = 34.9868366552d0 double precision, parameter :: dfour_pi = 4.d0*dacos(-1.d0) double precision, parameter :: dtwo_pi = 2.d0*dacos(-1.d0) +double precision, parameter :: inv_pi = 1.d0/dacos(-1.d0) double precision, parameter :: inv_sq_pi = 1.d0/dsqrt(dacos(-1.d0)) double precision, parameter :: inv_sq_pi_2 = 0.5d0/dsqrt(dacos(-1.d0)) double precision, parameter :: thresh = 1.d-15 From 227c139895efb89f12b0b84f2b93b71204babfae Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Tue, 2 Jun 2020 13:36:00 -0500 Subject: [PATCH 4/8] smaller three to four index transformation --- src/mo_two_e_ints/df_mo_ints.irp.f | 212 ++++++++++++++++++++++++ src/mo_two_e_ints/mo_bi_integrals.irp.f | 3 +- 2 files changed, 214 insertions(+), 1 deletion(-) diff --git a/src/mo_two_e_ints/df_mo_ints.irp.f b/src/mo_two_e_ints/df_mo_ints.irp.f index dbb10782..c9d03e0c 100644 --- a/src/mo_two_e_ints/df_mo_ints.irp.f +++ b/src/mo_two_e_ints/df_mo_ints.irp.f @@ -19,6 +19,218 @@ BEGIN_PROVIDER [complex*16, df_mo_integrals_complex, (mo_num_per_kpt,mo_num_per_ END_PROVIDER +subroutine mo_map_fill_from_df_single + use map_module + implicit none + BEGIN_DOC + ! fill mo bielec integral map using 3-index df integrals + END_DOC + + integer :: i,k,j,l + integer :: ki,kk,kj,kl + integer :: ii,ik,ij,il + integer :: kikk2,kjkl2,jl2,ik2 + integer :: i_mo,j_mo,i_df + + complex*16,allocatable :: ints_ik(:,:,:), ints_jl(:,:,:) + + complex*16 :: integral + integer :: n_integrals_1, n_integrals_2 + integer :: size_buffer + integer(key_kind),allocatable :: buffer_i_1(:), buffer_i_2(:) + real(integral_kind),allocatable :: buffer_values_1(:), buffer_values_2(:) + double precision :: tmp_re,tmp_im + integer :: mo_num_kpt_2 + + double precision :: cpu_1, cpu_2, wall_1, wall_2, wall_0 + double precision :: map_mb + + logical :: use_map1 + integer(key_kind) :: idx_tmp + double precision :: sign + + mo_num_kpt_2 = mo_num_per_kpt * mo_num_per_kpt + + size_buffer = min(mo_num_per_kpt*mo_num_per_kpt*mo_num_per_kpt,16000000) + print*, 'Providing the mo_bielec integrals from 3-index df integrals' + call write_time(6) +! call ezfio_set_integrals_bielec_disk_access_mo_integrals('Write') +! TOUCH read_mo_integrals read_ao_integrals write_mo_integrals write_ao_integrals + + call wall_time(wall_1) + call cpu_time(cpu_1) + + allocate( ints_jl(mo_num_per_kpt,mo_num_per_kpt,df_num)) + allocate( ints_ik(mo_num_per_kpt,mo_num_per_kpt,df_num)) + + wall_0 = wall_1 + do kl=1, kpt_num + do kj=1, kl + call idx2_tri_int(kj,kl,kjkl2) + if (kj < kl) then + do i_mo=1,mo_num_per_kpt + do j_mo=1,mo_num_per_kpt + do i_df=1,df_num + ints_jl(i_mo,j_mo,i_df) = dconjg(df_mo_integrals_complex(j_mo,i_mo,i_df,kjkl2)) + enddo + enddo + enddo + else + ints_jl = df_mo_integrals_complex(:,:,:,kjkl2) + endif + + 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_mo,j_mo,i_df) = dconjg(df_mo_integrals_complex(j_mo,i_mo,i_df,kikk2)) + enddo + enddo + enddo +! ints_ik = conjg(reshape(df_mo_integral_array(:,:,:,kikk2),(/mo_num_per_kpt,mo_num_per_kpt,df_num/),order=(/2,1,3/))) + else + ints_ik = df_mo_integrals_complex(:,:,:,kikk2) + endif + + !$OMP PARALLEL PRIVATE(i,k,j,l,ii,ik,ij,il,jl2,ik2, & + !$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 mu=1,df_num + 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) + mjl = ints_jl(ij,il,mu) + if (mjl.eq.(0.d0,0.d0)) cycle + 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 + mik = ints_ik(ii,ik,mu) + integral = mik * dconjg(mjl) +! print*,i,k,j,l,real(integral),imag(integral) + if (cdabs(integral) < mo_integrals_threshold) then + cycle + endif + call ao_two_e_integral_complex_map_idx_sign(i,j,k,l,use_map1,idx_tmp,sign) + tmp_re = dble(integral) + tmp_im = dimag(integral) + if (use_map1) then + n_integrals_1 += 1 + buffer_i_1(n_integrals_1)=idx_tmp + buffer_values_1(n_integrals_1)=tmp_re + if (sign.ne.0.d0) then + n_integrals_1 += 1 + buffer_i_1(n_integrals_1)=idx_tmp+1 + buffer_values_1(n_integrals_1)=tmp_im*sign + endif + if (n_integrals_1 >= size(buffer_i_1)-1) then + !call map_append(mo_integrals_map, buffer_i_1, buffer_values_1, n_integrals_1) + call insert_into_ao_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1) + n_integrals_1 = 0 + endif + else + n_integrals_2 += 1 + buffer_i_2(n_integrals_2)=idx_tmp + buffer_values_2(n_integrals_2)=tmp_re + if (sign.ne.0.d0) then + n_integrals_2 += 1 + buffer_i_2(n_integrals_2)=idx_tmp+1 + buffer_values_2(n_integrals_2)=tmp_im*sign + endif + if (n_integrals_2 >= size(buffer_i_2)-1) then + !call map_append(mo_integrals_map_2, buffer_i_2, buffer_values_2, n_integrals_2) + call insert_into_ao_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_2) + n_integrals_2 = 0 + endif + endif + + enddo !ii + enddo !ik + enddo !ij + enddo !il + enddo !mu + !$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_ao_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1) + endif + if (n_integrals_2 > 0) then + !call map_append(mo_integrals_map_2, buffer_i_2, buffer_values_2, n_integrals_2) + call insert_into_ao_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_2) + endif + 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_df_single + subroutine mo_map_fill_from_df use map_module implicit none 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 21422ba3..a4b3530e 100644 --- a/src/mo_two_e_ints/mo_bi_integrals.irp.f +++ b/src/mo_two_e_ints/mo_bi_integrals.irp.f @@ -43,7 +43,8 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ] return else if (read_df_mo_integrals.or.read_df_ao_integrals) then PROVIDE df_mo_integrals_complex - call mo_map_fill_from_df + !call mo_map_fill_from_df + call mo_map_fill_from_df_single return else PROVIDE ao_two_e_integrals_in_map From fc8abcbf0a081a7393edc79b9bcdcde92d9447b3 Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Tue, 2 Jun 2020 16:10:16 -0500 Subject: [PATCH 5/8] minor fix --- src/mo_two_e_ints/df_mo_ints.irp.f | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/mo_two_e_ints/df_mo_ints.irp.f b/src/mo_two_e_ints/df_mo_ints.irp.f index c9d03e0c..9bc1e3c0 100644 --- a/src/mo_two_e_ints/df_mo_ints.irp.f +++ b/src/mo_two_e_ints/df_mo_ints.irp.f @@ -26,7 +26,7 @@ subroutine mo_map_fill_from_df_single ! fill mo bielec integral map using 3-index df integrals END_DOC - integer :: i,k,j,l + integer :: i,k,j,l,mu integer :: ki,kk,kj,kl integer :: ii,ik,ij,il integer :: kikk2,kjkl2,jl2,ik2 @@ -34,7 +34,7 @@ subroutine mo_map_fill_from_df_single complex*16,allocatable :: ints_ik(:,:,:), ints_jl(:,:,:) - complex*16 :: integral + 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(:) @@ -97,6 +97,7 @@ subroutine mo_map_fill_from_df_single 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) & From 1c14b837c298f3122ac97214f45dea5b0a03a7e4 Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Tue, 2 Jun 2020 18:16:36 -0500 Subject: [PATCH 6/8] int type --- src/utils_complex/export_integrals_ao_cplx.irp.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils_complex/export_integrals_ao_cplx.irp.f b/src/utils_complex/export_integrals_ao_cplx.irp.f index d24a51e2..7d40ba76 100644 --- a/src/utils_complex/export_integrals_ao_cplx.irp.f +++ b/src/utils_complex/export_integrals_ao_cplx.irp.f @@ -159,7 +159,7 @@ provide ao_two_e_integrals_in_map if (cdabs(tmp_cmplx-int2e_tmp1).gt.thr0) then print'(4(I4),4(E15.7))',i,j,k,l,tmp_cmplx,int2e_tmp1 endif - integer*8 :: ii + integer :: ii ii = l-ao_integrals_cache_min ii = ior( shiftl(ii,6), k-ao_integrals_cache_min) ii = ior( shiftl(ii,6), j-ao_integrals_cache_min) From fee0ae8680fa36d376a69048e910acb9ebdd4c6c Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Tue, 2 Jun 2020 18:17:15 -0500 Subject: [PATCH 7/8] better mo transformation --- src/mo_two_e_ints/df_mo_ints.irp.f | 230 +++++++++++++++++++++++- src/mo_two_e_ints/mo_bi_integrals.irp.f | 3 +- 2 files changed, 228 insertions(+), 5 deletions(-) diff --git a/src/mo_two_e_ints/df_mo_ints.irp.f b/src/mo_two_e_ints/df_mo_ints.irp.f index 9bc1e3c0..3a61911e 100644 --- a/src/mo_two_e_ints/df_mo_ints.irp.f +++ b/src/mo_two_e_ints/df_mo_ints.irp.f @@ -19,6 +19,228 @@ BEGIN_PROVIDER [complex*16, df_mo_integrals_complex, (mo_num_per_kpt,mo_num_per_ END_PROVIDER +subroutine mo_map_fill_from_df_dot + use map_module + implicit none + BEGIN_DOC + ! fill mo bielec integral map using 3-index df 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 + + 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) +! 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_df_dot + subroutine mo_map_fill_from_df_single use map_module implicit none @@ -155,7 +377,7 @@ subroutine mo_map_fill_from_df_single endif if (n_integrals_1 >= size(buffer_i_1)-1) then !call map_append(mo_integrals_map, buffer_i_1, buffer_values_1, n_integrals_1) - call insert_into_ao_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1) + call insert_into_mo_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1,mo_integrals_threshold) n_integrals_1 = 0 endif else @@ -169,7 +391,7 @@ subroutine mo_map_fill_from_df_single endif if (n_integrals_2 >= size(buffer_i_2)-1) then !call map_append(mo_integrals_map_2, buffer_i_2, buffer_values_2, n_integrals_2) - call insert_into_ao_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_2) + 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 @@ -183,11 +405,11 @@ subroutine mo_map_fill_from_df_single if (n_integrals_1 > 0) then !call map_append(mo_integrals_map, buffer_i_1, buffer_values_1, n_integrals_1) - call insert_into_ao_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1) + 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_ao_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_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, & 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 a4b3530e..16322c19 100644 --- a/src/mo_two_e_ints/mo_bi_integrals.irp.f +++ b/src/mo_two_e_ints/mo_bi_integrals.irp.f @@ -44,7 +44,8 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ] else if (read_df_mo_integrals.or.read_df_ao_integrals) then PROVIDE df_mo_integrals_complex !call mo_map_fill_from_df - call mo_map_fill_from_df_single + !call mo_map_fill_from_df_single + call mo_map_fill_from_df_dot return else PROVIDE ao_two_e_integrals_in_map From 335386fa784bd04f1da38bf5c247395635e47891 Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Thu, 11 Jun 2020 13:32:24 -0500 Subject: [PATCH 8/8] fixed integral transformation; added complex fcidump; fixed kpts bitmasks --- src/bitmask/core_inact_act_virt.irp.f | 10 +-- src/mo_two_e_ints/df_mo_ints.irp.f | 6 +- src/tools/fcidump.irp.f | 91 +++++++++++++++++++++++++++ 3 files changed, 100 insertions(+), 7 deletions(-) diff --git a/src/bitmask/core_inact_act_virt.irp.f b/src/bitmask/core_inact_act_virt.irp.f index d2efef89..ae00d774 100644 --- a/src/bitmask/core_inact_act_virt.irp.f +++ b/src/bitmask/core_inact_act_virt.irp.f @@ -448,7 +448,7 @@ BEGIN_PROVIDER [ integer, n_core_orb_kpts, (kpt_num)] do k=1,kpt_num n_core_orb_kpts(k) = 0 - kshift = (1-k)*mo_num_per_kpt + kshift = (k-1)*mo_num_per_kpt do i = 1, mo_num_per_kpt if(mo_class(i+kshift) == 'Core')then n_core_orb_kpts(k) += 1 @@ -469,7 +469,7 @@ BEGIN_PROVIDER [ integer, n_inact_orb_kpts, (kpt_num)] do k=1,kpt_num n_inact_orb_kpts(k) = 0 - kshift = (1-k)*mo_num_per_kpt + kshift = (k-1)*mo_num_per_kpt do i = 1, mo_num_per_kpt if(mo_class(i+kshift) == 'Inactive')then n_inact_orb_kpts(k) += 1 @@ -490,7 +490,7 @@ BEGIN_PROVIDER [ integer, n_act_orb_kpts, (kpt_num)] do k=1,kpt_num n_act_orb_kpts(k) = 0 - kshift = (1-k)*mo_num_per_kpt + kshift = (k-1)*mo_num_per_kpt do i = 1, mo_num_per_kpt if(mo_class(i+kshift) == 'Active')then n_act_orb_kpts(k) += 1 @@ -511,7 +511,7 @@ BEGIN_PROVIDER [ integer, n_virt_orb_kpts, (kpt_num)] do k=1,kpt_num n_virt_orb_kpts(k) = 0 - kshift = (1-k)*mo_num_per_kpt + kshift = (k-1)*mo_num_per_kpt do i = 1, mo_num_per_kpt if(mo_class(i+kshift) == 'Virtual')then n_virt_orb_kpts(k) += 1 @@ -532,7 +532,7 @@ BEGIN_PROVIDER [ integer, n_del_orb_kpts, (kpt_num)] do k=1,kpt_num n_del_orb_kpts(k) = 0 - kshift = (1-k)*mo_num_per_kpt + kshift = (k-1)*mo_num_per_kpt do i = 1, mo_num_per_kpt if(mo_class(i+kshift) == 'Deleted')then n_del_orb_kpts(k) += 1 diff --git a/src/mo_two_e_ints/df_mo_ints.irp.f b/src/mo_two_e_ints/df_mo_ints.irp.f index 3a61911e..eba3b3da 100644 --- a/src/mo_two_e_ints/df_mo_ints.irp.f +++ b/src/mo_two_e_ints/df_mo_ints.irp.f @@ -48,7 +48,8 @@ subroutine mo_map_fill_from_df_dot logical :: use_map1 integer(key_kind) :: idx_tmp double precision :: sign - complex*16, external :: zdotc + !complex*16, external :: zdotc + complex*16, external :: zdotu mo_num_kpt_2 = mo_num_per_kpt * mo_num_per_kpt @@ -145,7 +146,8 @@ subroutine mo_map_fill_from_df_dot 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 = 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 diff --git a/src/tools/fcidump.irp.f b/src/tools/fcidump.irp.f index bf4d07fb..de878dc6 100644 --- a/src/tools/fcidump.irp.f +++ b/src/tools/fcidump.irp.f @@ -18,6 +18,97 @@ program fcidump ! electrons ! END_DOC + if (is_complex) then + call fcidump_complex + else + call fcidump_real + endif +end + +subroutine fcidump_complex + implicit none + character*(128) :: output + integer :: i_unit_output,getUnitAndOpen + output=trim(ezfio_filename)//'.FCIDUMP' + i_unit_output = getUnitAndOpen(output,'w') + + integer :: i,j,k,l + integer :: i1,j1,k1,l1 + integer :: i2,j2,k2,l2,ik2,jl2 + integer :: ki,kj,kk,kl + integer :: ii,ij,ik,il + integer*8 :: m + character*(2), allocatable :: A(:) + + write(i_unit_output,*) '&FCI NORB=', n_act_orb, ', NELEC=', elec_num-n_core_orb*2, & + ', MS2=', (elec_alpha_num-elec_beta_num), ',' + allocate (A(n_act_orb)) + A = '1,' + write(i_unit_output,*) 'ORBSYM=', (A(i), i=1,n_act_orb) + write(i_unit_output,*) 'ISYM=0,' + write(i_unit_output,*) '/' + deallocate(A) + + integer(key_kind), allocatable :: keys(:) + double precision, allocatable :: values(:) + integer(cache_map_size_kind) :: n_elements, n_elements_max + PROVIDE mo_two_e_integrals_in_map + + complex*16 :: get_two_e_integral_complex, integral + + do kl=1,kpt_num + do kj=1,kl + do kk=1,kl + ki=kconserv(kl,kk,kj) + if (ki>kl) cycle + do l1=1,n_act_orb_kpts(kl) + il=list_act_kpts(l1,kl) + l = (kl-1)*mo_num_per_kpt + il + do j1=1,n_act_orb_kpts(kj) + ij=list_act_kpts(j1,kj) + j = (kj-1)*mo_num_per_kpt + ij + if (j>l) exit + call idx2_tri_int(j,l,jl2) + do k1=1,n_act_orb_kpts(kk) + ik=list_act_kpts(k1,kk) + k = (kk-1)*mo_num_per_kpt + ik + if (k>l) exit + do i1=1,n_act_orb_kpts(ki) + ii=list_act_kpts(i1,ki) + i = (ki-1)*mo_num_per_kpt + ii + if ((j==l) .and. (i>k)) exit + call idx2_tri_int(i,k,ik2) + if (ik2 > jl2) exit + integral = get_two_e_integral_complex(i,j,k,l,mo_integrals_map,mo_integrals_map_2) + if (cdabs(integral) > mo_integrals_threshold) then + write(i_unit_output,'(2(E25.15,X),4(I6,X))') dble(integral), dimag(integral),i,k,j,l + endif + enddo + enddo + enddo + enddo + enddo + enddo + enddo + + do kj=1,kpt_num + do j1=1,n_act_orb_kpts(kj) + ij = list_act_kpts(j1,kj) + j = (kj-1)*mo_num_per_kpt + ij + do i1=j1,n_act_orb_kpts(kj) + ii = list_act_kpts(i1,kj) + i = (kj-1)*mo_num_per_kpt + ii + integral = mo_one_e_integrals_kpts(ii,ij,kj) + core_fock_operator_complex(i,j) + if (cdabs(integral) > mo_integrals_threshold) then + write(i_unit_output,'(2(E25.15,X),4(I6,X))') dble(integral),dimag(integral), i,j,0,0 + endif + enddo + enddo + enddo + write(i_unit_output,*) core_energy, 0, 0, 0, 0 +end +subroutine fcidump_real + implicit none character*(128) :: output integer :: i_unit_output,getUnitAndOpen output=trim(ezfio_filename)//'.FCIDUMP'