program psiSVD_naiv1by1_v1 implicit none BEGIN_DOC ! perturbative approach to build psi_postsvd END_DOC read_wf = .True. TOUCH read_wf PROVIDE N_int call run() end subroutine run USE OMP_LIB implicit none integer(bit_kind) :: det1(N_int,2), det2(N_int,2) integer :: degree, i_state integer :: i, j, k, l, m, n double precision :: x, y, h12 double precision, allocatable :: Uref(:,:), Dref(:), Vtref(:,:), Aref(:,:), Vref(:,:) integer :: rank_max double precision :: E0, overlop, Ept2, Em, norm double precision, allocatable :: H0(:,:) double precision, allocatable :: eigvec0(:,:), eigval0(:), coeff_psi(:), coeff_tmp(:) integer :: ii, jj, ia, ib, ja, jb double precision, allocatable :: Hdiag(:), H0_1d(:), H0_tmp(:,:) integer :: na_new, nb_new, ind_new, ind_gs double precision :: ctmp, coeff_new double precision, allocatable :: epsil(:), epsil_energ(:), check_ov(:) double precision, allocatable :: Uezfio(:,:,:), Dezfio(:,:), Vezfio(:,:,:) integer :: n_selected, n_toselect, n_tmp, na_max, nb_max integer, allocatable :: numalpha_selected(:), numbeta_selected(:) integer, allocatable :: numalpha_tmp(:), numbeta_tmp(:) double precision :: t_beg, t_end integer(kind=8) :: W_tbeg, W_tend, W_tbeg_it, W_tend_it, W_ir real(kind=8) :: W_tot_time, W_tot_time_it integer :: nb_taches !$OMP PARALLEL nb_taches = OMP_GET_NUM_THREADS() !$OMP END PARALLEL call SYSTEM_CLOCK(COUNT=W_tbeg, COUNT_RATE=W_ir) i_state = 1 det1(:,1) = psi_det_alpha_unique(:,1) det2(:,1) = psi_det_alpha_unique(:,1) det1(:,2) = psi_det_beta_unique(:,1) det2(:,2) = psi_det_beta_unique(:,1) call get_excitation_degree_spin(det1(1,1),det2(1,1),degree,N_int) call get_excitation_degree(det1,det2,degree,N_int) call i_H_j(det1, det2, N_int, h12) ! --------------------------------------------------------------------------------------- ! construct the initial CISD matrix print *, ' ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~' print *, ' CISD matrix:', n_det_alpha_unique,'x',n_det_beta_unique print *, ' N det :', N_det print *, ' ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~' allocate( Aref(n_det_alpha_unique,n_det_beta_unique) ) Aref(:,:) = 0.d0 do k = 1, N_det i = psi_bilinear_matrix_rows(k) j = psi_bilinear_matrix_columns(k) Aref(i,j) = psi_bilinear_matrix_values(k,i_state) enddo ! --------------------------------------------------------------------------------------- ! --------------------------------------------------------------------------------------- ! perform a Full SVD allocate( Uref(n_det_alpha_unique,n_det_alpha_unique) ) allocate( Dref(min(n_det_alpha_unique,n_det_beta_unique)) ) allocate( Vtref(n_det_beta_unique,n_det_beta_unique) ) call cpu_time(t_beg) call svd_s(Aref, size(Aref,1), Uref, size(Uref,1), Dref, Vtref & , size(Vtref,1), n_det_alpha_unique, n_det_beta_unique) call cpu_time(t_end) print *, " SVD is performed after (min)", (t_end-t_beg)/60. allocate( Vref(n_det_beta_unique,n_det_beta_unique) ) do l = 1, n_det_beta_unique do i = 1, n_det_beta_unique Vref(i,l) = Vtref(l,i) enddo enddo deallocate( Vtref ) deallocate( Aref ) ! --------------------------------------------------------------------------------------- ! *** PARAMETERS *** ! na_max = n_det_alpha_unique nb_max = n_det_beta_unique ! *** ***** *** ! print *, ' na_max = ', na_max print *, ' nb_max = ', nb_max ! --------------------------------------------------------------------------------------- ! initial wavefunction: psi_0 n_selected = 1 allocate( numalpha_selected(n_selected) , numbeta_selected(n_selected) ) numalpha_selected(1) = 1 numbeta_selected (1) = 1 ! --------------------------------------------------------------------------------------- ! --------------------------------------------------------------------------------------- ! construnc the initial basis to select phi_1 from the FSVD n_toselect = min(na_max,nb_max) - n_selected print *, ' toselect = ', n_toselect print *, ' to trun = ', n_det_alpha_unique*n_det_beta_unique - na_max*nb_max ! --------------------------------------------------------------------------------------- ! --------------------------------------------------------------------------------------- ! read < u_k v_l | H | u_k v_l > for all vectors allocate( Hdiag(n_toselect) , H0(n_selected,n_selected) ) n_tmp = n_det_alpha_unique * n_det_beta_unique - 1 open( unit=11, FILE="klHkl_v1.dat", ACTION="READ") !open( unit=11, FILE="klHkl_v2.dat", ACTION="READ") read(11,*) i, i, E0 H0(1,1) = E0 do i = 1, n_tmp read(11,*) ia, ib, ctmp if( ia .eq. ib ) then ii = ia - 1 Hdiag(ii) = ctmp !print *, ia, ib , Hdiag(ia-1) endif enddo close(11) ! --------------------------------------------------------------------------------------- E0 = E0 + nuclear_repulsion Em = E0 print*, ' space dimen = ', n_selected print*, ' ground state Em = ', Em print*, ' ground state E0 = ', E0 na_new = 1 nb_new = 1 ind_new = 0 !________________________________________________________________________________________________________ ! ! increase the size of psi0 iteratively !________________________________________________________________________________________________________ ! *** PARAMETERS *** ! rank_max = min( na_max , nb_max ) - 1 ! *** ***** *** ! if( rank_max .gt. n_toselect ) then print *, " rank_max should be less then n_toselect" stop endif do while( n_selected .lt. rank_max ) call SYSTEM_CLOCK(COUNT=W_tbeg_it, COUNT_RATE=W_ir) print*, ' ' print*, ' new iteration ' if( n_toselect .lt. 1 ) then print*, ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ' print*, ' no more vectors to construct a new basis ' print*, ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ' stop else ! --------------------------------------------------------------------------------------- ! select a new vector na_new += 1 nb_new += 1 ind_new += 1 print *, ' best vector', na_new, nb_new ! < psi_old | H | delta_psi > allocate( H0_1d(n_selected) ) call const_H0_1d(na_new, nb_new, na_max, nb_max, n_selected, Uref, Vref, numalpha_selected, numbeta_selected, H0_1d) ! --------------------------------------------------------------------------------------- ! --------------------------------------------------------------------------------------- ! new psi allocate( numalpha_tmp(n_selected), numbeta_tmp(n_selected) ) allocate( H0_tmp(n_selected,n_selected) ) numalpha_tmp(:) = numalpha_selected(:) numbeta_tmp (:) = numbeta_selected (:) H0_tmp (:,:) = H0 (:,:) deallocate( numalpha_selected, numbeta_selected, H0 ) n_tmp = n_selected n_selected = n_selected + 1 allocate( numalpha_selected(n_selected) , numbeta_selected(n_selected) ) allocate( H0(n_selected,n_selected) ) H0(:,:) = 0.d0 do l = 1, n_tmp numalpha_selected(l) = numalpha_tmp(l) numbeta_selected (l) = numbeta_tmp (l) enddo H0(1:n_tmp,1:n_tmp) = H0_tmp(1:n_tmp,1:n_tmp) deallocate( numalpha_tmp, numbeta_tmp, H0_tmp ) numalpha_selected(n_selected) = na_new numbeta_selected (n_selected) = nb_new H0(1:n_tmp,n_selected) = H0_1d(1:n_tmp) H0(n_selected,1:n_tmp) = H0_1d(1:n_tmp) deallocate( H0_1d ) H0(n_selected,n_selected) = Hdiag(ind_new) ! --------------------------------------------------------------------------------------- ! --------------------------------------------------------------------------------------- ! energy without diag ! < psi | psi > norm = 0.d0 do j = 1, n_selected ja = numalpha_selected(j) jb = numbeta_selected (j) if(ja.eq.jb) norm = norm + Dref(ja)*Dref(jb) enddo ! < psi | H | psi > Em = 0.d0 do j = 1, n_selected ja = numalpha_selected(j) jb = numbeta_selected (j) if(ja.eq.jb) then do i = 1, n_selected ia = numalpha_selected(i) ib = numbeta_selected (i) if(ia.eq.ib) Em = Em + Dref(ja) * H0(j,i) * Dref(ia) enddo endif enddo ! Em = < psi | H | psi > / < psi | psi > Em = Em / norm + nuclear_repulsion ! --------------------------------------------------------------------------------------- ! --------------------------------------------------------------------------------------- ! energy with diag allocate( eigvec0(n_selected,n_selected), eigval0(n_selected) ) call lapack_diag(eigval0, eigvec0, H0, n_selected, n_selected) ! get the postsvd ground state allocate( check_ov(n_selected) ) do l = 1, n_selected overlop = 0.d0 do i = 1, n_selected ia = numalpha_selected(i) ib = numbeta_selected (i) if( ia .eq. ib ) overlop = overlop + eigvec0(i,l) * Dref(ia) enddo check_ov(l) = dabs(overlop) enddo ind_gs = MAXLOC( check_ov, DIM=1 ) overlop = check_ov(ind_gs) E0 = eigval0(ind_gs)+nuclear_repulsion coeff_psi(:) = eigvec0(:,ind_gs) deallocate( check_ov, eigval0, eigvec0 ) ! --------------------------------------------------------------------------------------- print*, ' space dimen = ', n_selected print*, ' E bef diag = ', Em print*, ' E aft diag = ', E0 print*, ' overlop = ', overlop print*, ' index = ', ind_gs write(211, '( 3(I5,3X), 4(F15.8,3X) )') n_selected, na_new, nb_new, Em, E0, overlop ! --------------------------------------------------------------------------------------- ! remove selected pair | na_new nb_new > n_toselect = n_toselect - 1 print*, ' rank to select = ', n_toselect ! --------------------------------------------------------------------------------------- endif call SYSTEM_CLOCK(COUNT=W_tend_it, COUNT_RATE=W_ir) W_tot_time_it = real(W_tend_it-W_tbeg_it, kind=8) / real(W_ir, kind=8) print*, " " print*, " elapsed time (min) = ", W_tot_time_it/60.d0 end do !________________________________________________________________________________________________________ !________________________________________________________________________________________________________ ! *************************************************************************************************** ! save to ezfion !allocate( Uezfio(n_det_alpha_unique,rank0,1), Dezfio(rank0,1), Vezfio(n_det_beta_unique,rank0,1) ) !do l = 1, rank0 ! Dezfio(l,1) = coeff_psi(l) ! Uezfio(:,l,1) = U0(:,l) ! Vezfio(:,l,1) = V0(:,l) !enddo !call ezfio_set_spindeterminants_n_det(N_det) !call ezfio_set_spindeterminants_n_states(N_states) !call ezfio_set_spindeterminants_n_det_alpha(n_det_alpha_unique) !call ezfio_set_spindeterminants_n_det_beta(n_det_beta_unique) !call ezfio_set_spindeterminants_psi_coef_matrix_rows(psi_bilinear_matrix_rows) !call ezfio_set_spindeterminants_psi_coef_matrix_columns(psi_bilinear_matrix_columns) !call ezfio_set_spindeterminants_psi_coef_matrix_values(psi_bilinear_matrix_values) !call ezfio_set_spindeterminants_n_svd_coefs(rank0) !call ezfio_set_spindeterminants_psi_svd_alpha(Uezfio) !call ezfio_set_spindeterminants_psi_svd_beta(Vezfio ) !call ezfio_set_spindeterminants_psi_svd_coefs(Dezfio) !deallocate( Uezfio, Dezfio, Vezfio ) ! *************************************************************************************************** call SYSTEM_CLOCK(COUNT=W_tend, COUNT_RATE=W_ir) W_tot_time = real(W_tend - W_tbeg, kind=8) / real(W_ir, kind=8) print *, ' ___________________________________________________________________' print *, ' ' print *, " Execution avec ", nb_taches, " threads" print *, " total elapsed time (min) = ", W_tot_time/60.d0 print *, ' ___________________________________________________________________' deallocate( Dref ) deallocate( Uref, Vref ) deallocate( coeff_psi ) deallocate( numalpha_selected, numbeta_selected ) deallocate( H0, Hdiag ) end subroutine const_H0_1d(na_new, nb_new, na_max, nb_max, n_selected, Uref, Vref, numalpha_selected, numbeta_selected, H0_1d) implicit none integer, intent(in) :: na_new, nb_new, na_max, nb_max, n_selected integer, intent(in) :: numalpha_selected(n_selected), numbeta_selected(n_selected) double precision, intent(in) :: Uref(n_det_alpha_unique,n_det_alpha_unique) double precision, intent(in) :: Vref(n_det_beta_unique ,n_det_beta_unique) double precision, intent(out) :: H0_1d(n_selected) integer(bit_kind) :: det1(N_int,2) integer(bit_kind) :: det2(N_int,2) integer :: degree, na, nb integer :: i, j, k, l, ii, jj, m double precision :: h12 double precision, allocatable :: Hmat_kl(:,:), tmp1(:,:), tmp2(:,:) double precision, allocatable :: U1d(:), V1d(:) double precision, allocatable :: Utmp(:,:), Vtmp(:,:) double precision :: ti, tf print *, "" print *, " start const_H0_1d" call wall_time(ti) na = n_det_alpha_unique nb = n_det_beta_unique allocate( U1d(na) , V1d(nb) ) U1d(1:na) = Uref(1:na,na_new) V1d(1:nb) = Vref(1:nb,nb_new) allocate( tmp1(na,nb) ) tmp1 = 0.d0 !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(i,j,k,l,h12,det1,det2,degree,tmp2) & !$OMP SHARED(na,nb,psi_det_alpha_unique,psi_det_beta_unique, & !$OMP N_int,U1d,V1d,tmp1) allocate( tmp2(na,nb) ) tmp2 = 0.d0 !$OMP DO do l = 1, nb det2(:,2) = psi_det_beta_unique(:,l) do j = 1, nb det1(:,2) = psi_det_beta_unique(:,j) call get_excitation_degree_spin(det1(1,2),det2(1,2),degree,N_int) if(degree .gt. 2) cycle do k = 1, na det2(:,1) = psi_det_alpha_unique(:,k) do i = 1, na det1(:,1) = psi_det_alpha_unique(:,i) call get_excitation_degree(det1,det2,degree,N_int) if(degree .gt. 2) cycle call i_H_j(det1, det2, N_int, h12) if( h12 .eq. 0.d0) cycle tmp2(i,j) += h12 * U1d(k) * V1d(l) enddo enddo enddo enddo !$OMP END DO !$OMP CRITICAL do j = 1, nb do i = 1, na tmp1(i,j) += tmp2(i,j) enddo enddo !$OMP END CRITICAL deallocate( tmp2 ) !$OMP END PARALLEL deallocate( U1d , V1d ) ! tmp2(j,m) = sum_i tmp1(i,j) x Uref(i,m) allocate( Utmp(na,na_max) ) Utmp(1:na,1:na_max) = Uref(1:na,1:na_max) allocate( tmp2(nb,na_max) ) call DGEMM('T', 'N', nb, na_max, na, 1.d0, & tmp1, size(tmp1,1), Utmp, size(Utmp,1), & 0.d0, tmp2, size(tmp2,1) ) deallocate( tmp1 ) deallocate( Utmp ) ! Hmat_kl(m,n) = sum_j tmp2(j,m) x Vref(j,n) allocate( Vtmp(nb,nb_max) ) Vtmp(1:nb,1:nb_max) = Vref(1:nb,1:nb_max) allocate( Hmat_kl(na_max,nb_max) ) call DGEMM('T', 'N', na_max, nb_max, nb, 1.d0, & tmp2, size(tmp2,1), Vtmp, size(Vtmp,1), & 0.d0, Hmat_kl, size(Hmat_kl,1) ) deallocate( tmp2 ) deallocate( Vtmp ) do m = 1, n_selected ii = numalpha_selected(m) jj = numbeta_selected (m) H0_1d(m) = Hmat_kl(ii,jj) enddo deallocate( Hmat_kl ) call wall_time(tf) print *, " end const_H0_1d after (min) ", (tf-ti)/60. print *, "" return end subroutine const_H0_1d