subroutine CISD_SC2(dets_in,u_in,energies,diag_H_elements,dim_in,sze,N_st,Nint,convergence) use bitmasks implicit none BEGIN_DOC ! CISD+SC2 method :: take off all the disconnected terms of a CISD (selected or not) ! ! dets_in : bitmasks corresponding to determinants ! ! u_in : guess coefficients on the various states. Overwritten ! on exit ! ! dim_in : leftmost dimension of u_in ! ! sze : Number of determinants ! ! N_st : Number of eigenstates ! ! Initial guess vectors are not necessarily orthonormal END_DOC integer, intent(in) :: dim_in, sze, N_st, Nint integer(bit_kind), intent(in) :: dets_in(Nint,2,sze) double precision, intent(inout) :: u_in(dim_in,N_st) double precision, intent(out) :: energies(N_st) double precision, intent(out) :: diag_H_elements(dim_in) double precision, intent(in) :: convergence ASSERT (N_st > 0) ASSERT (sze > 0) ASSERT (Nint > 0) ASSERT (Nint == N_int) integer :: iter integer :: i,j,k,l,m logical :: converged double precision :: overlap(N_st,N_st) double precision :: u_dot_v, u_dot_u integer :: degree,N_double,index_hf double precision :: hij_elec, e_corr_double,e_corr,diag_h_mat_elem,inv_c0 double precision :: e_corr_double_before,accu,cpu_2,cpu_1 integer,allocatable :: degree_exc(:), index_double(:) integer :: i_ok double precision,allocatable :: e_corr_array(:),H_jj_ref(:),H_jj_dressed(:),hij_double(:) integer(bit_kind), allocatable :: doubles(:,:,:) allocate (doubles(Nint,2,sze),e_corr_array(sze),H_jj_ref(sze),H_jj_dressed(sze),& index_double(sze), degree_exc(sze), hij_double(sze)) call write_time(6) write(6,'(A)') '' write(6,'(A)') 'CISD SC2' write(6,'(A)') '========' !$OMP PARALLEL DEFAULT(NONE) & !$OMP SHARED(sze,N_st, & !$OMP H_jj_ref,Nint,dets_in,u_in) & !$OMP PRIVATE(i) !$OMP DO SCHEDULE(guided) do i=1,sze H_jj_ref(i) = diag_h_mat_elem(dets_in(1,1,i),Nint) enddo !$OMP END DO NOWAIT !$OMP END PARALLEL N_double = 0 e_corr = 0.d0 e_corr_double = 0.d0 do i = 1, sze call get_excitation_degree(ref_bitmask,dets_in(1,1,i),degree,Nint) degree_exc(i) = degree+1 if(degree==0)then index_hf=i else if (degree == 2)then N_double += 1 index_double(N_double) = i doubles(:,:,N_double) = dets_in(:,:,i) call i_H_j(ref_bitmask,dets_in(1,1,i),Nint,hij_elec) hij_double(N_double) = hij_elec e_corr_array(N_double) = u_in(i,1)* hij_elec e_corr_double += e_corr_array(N_double) e_corr += e_corr_array(N_double) else if (degree == 1)then call i_H_j(ref_bitmask,dets_in(1,1,i),Nint,hij_elec) e_corr += u_in(i,1)* hij_elec endif enddo inv_c0 = 1.d0/u_in(index_hf,1) do i = 1, N_double e_corr_array(i) = e_corr_array(i) * inv_c0 enddo e_corr = e_corr * inv_c0 e_corr_double = e_corr_double * inv_c0 converged = .False. e_corr_double_before = e_corr_double iter = 0 do while (.not.converged) iter +=1 !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(i,j,degree,accu) & !$OMP SHARED(H_jj_dressed,sze,H_jj_ref,index_hf,N_int,N_double,& !$OMP dets_in,doubles,degree_exc,e_corr_array,e_corr_double) !$OMP DO SCHEDULE(STATIC) do i=1,sze H_jj_dressed(i) = H_jj_ref(i) if (i==index_hf)cycle accu = -e_corr_double select case (N_int) case (1) do j=1,N_double degree = & popcnt(xor( dets_in(1,1,i),doubles(1,1,j))) + & popcnt(xor( dets_in(1,2,i),doubles(1,2,j))) if (degree<=ishft(degree_exc(i),1)) then accu += e_corr_array(j) endif enddo case (2) do j=1,N_double degree = & popcnt(xor( dets_in(1,1,i),doubles(1,1,j))) + & popcnt(xor( dets_in(1,2,i),doubles(1,2,j))) + & popcnt(xor( dets_in(2,1,i),doubles(2,1,j))) + & popcnt(xor( dets_in(2,2,i),doubles(2,2,j))) if (degree<=ishft(degree_exc(i),1)) then accu += e_corr_array(j) endif enddo case (3) do j=1,N_double degree = & popcnt(xor( dets_in(1,1,i),doubles(1,1,j))) + & popcnt(xor( dets_in(1,2,i),doubles(1,2,j))) + & popcnt(xor( dets_in(2,1,i),doubles(2,1,j))) + & popcnt(xor( dets_in(2,2,i),doubles(2,2,j))) + & popcnt(xor( dets_in(3,1,i),doubles(3,1,j))) + & popcnt(xor( dets_in(3,2,i),doubles(3,2,j))) if (degree<=ishft(degree_exc(i),1)) then accu += e_corr_array(j) endif enddo case default do j=1,N_double call get_excitation_degree(dets_in(1,1,i),doubles(1,1,j),degree,N_int) if (degree<=degree_exc(i)) then accu += e_corr_array(j) endif enddo end select H_jj_dressed(i) -= accu enddo !$OMP END DO !$OMP END PARALLEL if(sze<=N_det_max_jacobi)then double precision, allocatable :: eigenvectors(:,:), eigenvalues(:),H_matrix_tmp(:,:) allocate (H_matrix_tmp(size(H_matrix_all_dets,1),sze),eigenvalues(sze),eigenvectors(size(H_matrix_all_dets,1),sze)) do j=1,sze do i=1,sze H_matrix_tmp(i,j) = H_matrix_all_dets(i,j) enddo enddo do i = 1,sze H_matrix_tmp(i,i) = H_jj_dressed(i) enddo call lapack_diag(eigenvalues,eigenvectors, & H_matrix_tmp,size(H_matrix_all_dets,1),sze) do j=1,min(N_states_diag,sze) do i=1,sze u_in(i,j) = eigenvectors(i,j) enddo energies(j) = eigenvalues(j) enddo deallocate (H_matrix_tmp, eigenvalues, eigenvectors) else call davidson_diag_hjj(dets_in,u_in,H_jj_dressed,energies,dim_in,sze,N_st,N_states_diag,Nint,6) endif e_corr_double = 0.d0 inv_c0 = 1.d0/u_in(index_hf,1) do i = 1, N_double e_corr_array(i) = u_in(index_double(i),1)*inv_c0 * hij_double(i) e_corr_double += e_corr_array(i) enddo write(6,'(A,I3)') 'SC2 Iteration ', iter write(6,'(A)') '------------------' write(6,'(A)') '' write(6,'(A)') '===== ================' write(6,'(A)') 'State Energy ' write(6,'(A)') '===== ================' do i=1,N_st write(6,'(I5,1X,F16.10)') i, energies(i)+nuclear_repulsion enddo write(6,'(A)') '===== ================' write(6,'(A)') '' call write_double(6,(e_corr_double - e_corr_double_before),& 'Delta(E_corr)') converged = dabs(e_corr_double - e_corr_double_before) < convergence converged = converged if (converged) then do i = 1, dim_in diag_H_elements(i) = H_jj_dressed(i) - H_jj_ref(i) enddo exit endif e_corr_double_before = e_corr_double enddo call write_time(6) deallocate (doubles,e_corr_array,H_jj_ref,H_jj_dressed, & index_double, degree_exc, hij_double) end