diff --git a/plugins/FOBOCI/routines_foboci.irp.f b/plugins/FOBOCI/routines_foboci.irp.f index 4aca60d7..6fc60fae 100644 --- a/plugins/FOBOCI/routines_foboci.irp.f +++ b/plugins/FOBOCI/routines_foboci.irp.f @@ -799,7 +799,7 @@ end call dress_diag_elem_2h2p(dressing_H_mat_elem,N_det) call dress_diag_elem_2h1p(dressing_H_mat_elem,N_det,lmct,i_hole) call dress_diag_elem_1h2p(dressing_H_mat_elem,N_det,lmct,i_hole) - call davidson_diag_hjj(psi_det,psi_coef,dressing_H_mat_elem,energies,size(psi_coef,1),N_det,N_states_diag,N_int,output_determinants) + call davidson_diag_hjj(psi_det,psi_coef,dressing_H_mat_elem,energies,size(psi_coef,1),N_det,N_states,N_states_diag,N_int,output_determinants) do i = 1, 2 print*,'psi_coef = ',psi_coef(i,1) enddo diff --git a/plugins/MRCC_Utils/davidson.irp.f b/plugins/MRCC_Utils/davidson.irp.f index 713b3f61..7bde219a 100644 --- a/plugins/MRCC_Utils/davidson.irp.f +++ b/plugins/MRCC_Utils/davidson.irp.f @@ -1,7 +1,4 @@ - - - -subroutine davidson_diag_mrcc(dets_in,u_in,energies,dim_in,sze,N_st,Nint,iunit,istate) +subroutine davidson_diag_mrcc(dets_in,u_in,energies,dim_in,sze,N_st,N_st_diag,Nint,iunit,istate) use bitmasks implicit none BEGIN_DOC @@ -22,7 +19,7 @@ subroutine davidson_diag_mrcc(dets_in,u_in,energies,dim_in,sze,N_st,Nint,iunit,i ! ! Initial guess vectors are not necessarily orthonormal END_DOC - integer, intent(in) :: dim_in, sze, N_st, Nint, iunit, istate + integer, intent(in) :: dim_in, sze, N_st, Nint, iunit, istate, N_st_diag 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) @@ -31,6 +28,7 @@ subroutine davidson_diag_mrcc(dets_in,u_in,energies,dim_in,sze,N_st,Nint,iunit,i double precision :: diag_h_mat_elem integer :: i ASSERT (N_st > 0) + ASSERT (N_st_diag >= N_st) ASSERT (sze > 0) ASSERT (Nint > 0) ASSERT (Nint == N_int) @@ -52,11 +50,11 @@ subroutine davidson_diag_mrcc(dets_in,u_in,energies,dim_in,sze,N_st,Nint,iunit,i !$OMP END DO !$OMP END PARALLEL - call davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iunit,istate) + call davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag,Nint,iunit,istate) deallocate (H_jj) end -subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iunit,istate) +subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag,Nint,iunit,istate) use bitmasks implicit none BEGIN_DOC @@ -79,7 +77,7 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nin ! ! Initial guess vectors are not necessarily orthonormal END_DOC - integer, intent(in) :: dim_in, sze, N_st, Nint, istate + integer, intent(in) :: dim_in, sze, N_st, Nint, istate, N_st_diag integer(bit_kind), intent(in) :: dets_in(Nint,2,sze) double precision, intent(in) :: H_jj(sze) integer, intent(in) :: iunit @@ -106,6 +104,9 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nin double precision :: cpu, wall !PROVIDE det_connections +if (N_st_diag /= N_st) then + stop 'N_st_diag /= N_st todo in davidson' +endif call write_time(iunit) call wall_time(wall) diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index 91e8c120..9876c756 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -149,7 +149,7 @@ END_PROVIDER if (diag_algorithm == "Davidson") then call davidson_diag_mrcc(psi_det,CI_eigenvectors_dressed,CI_electronic_energy_dressed,& - size(CI_eigenvectors_dressed,1),N_det,N_states_diag,N_int,output_determinants,mrcc_state) + size(CI_eigenvectors_dressed,1),N_det,N_states,N_states_diag,N_int,output_determinants,mrcc_state) call u_0_S2_u_0(CI_eigenvectors_s2_dressed,CI_eigenvectors_dressed,N_det,psi_det,N_int,& N_states_diag,size(CI_eigenvectors_dressed,1)) diff --git a/plugins/MRCC_Utils_new/EZFIO.cfg b/plugins/MRCC_Utils_new/EZFIO.cfg deleted file mode 100644 index 789f30ef..00000000 --- a/plugins/MRCC_Utils_new/EZFIO.cfg +++ /dev/null @@ -1,4 +0,0 @@ -[energy] -type: double precision -doc: Calculated MRCC energy -interface: ezfio \ No newline at end of file diff --git a/plugins/MRCC_Utils_new/NEEDED_CHILDREN_MODULES b/plugins/MRCC_Utils_new/NEEDED_CHILDREN_MODULES deleted file mode 100644 index 5b16423e..00000000 --- a/plugins/MRCC_Utils_new/NEEDED_CHILDREN_MODULES +++ /dev/null @@ -1 +0,0 @@ -Perturbation Selectors_full Generators_full Psiref_Utils diff --git a/plugins/MRCC_Utils_new/README.rst b/plugins/MRCC_Utils_new/README.rst deleted file mode 100644 index 6f070867..00000000 --- a/plugins/MRCC_Utils_new/README.rst +++ /dev/null @@ -1,168 +0,0 @@ -=========== -MRCC Module -=========== - -Multi-Reference Coupled Cluster. - -Needed Modules -============== - -.. Do not edit this section. It was auto-generated from the -.. by the `update_README.py` script. - -.. image:: tree_dependency.png - -* `Perturbation `_ -* `Selectors_full `_ -* `Generators_full `_ -* `Psiref_Utils `_ - -Documentation -============= - -.. Do not edit this section. It was auto-generated from the -.. by the `update_README.py` script. - -`apply_excitation_operator `_ - Undocumented - - -`ci_eigenvectors_dressed `_ - Eigenvectors/values of the CI matrix - - -`ci_eigenvectors_s2_dressed `_ - Eigenvectors/values of the CI matrix - - -`ci_electronic_energy_dressed `_ - Eigenvectors/values of the CI matrix - - -`ci_energy_dressed `_ - N_states lowest eigenvalues of the dressed CI matrix - - -`davidson_diag_hjj_mrcc `_ - Davidson diagonalization with specific diagonal elements of the H matrix - .br - H_jj : specific diagonal H matrix elements to diagonalize de Davidson - .br - dets_in : bitmasks corresponding to determinants - .br - u_in : guess coefficients on the various states. Overwritten - on exit - .br - dim_in : leftmost dimension of u_in - .br - sze : Number of determinants - .br - N_st : Number of eigenstates - .br - iunit : Unit for the I/O - .br - Initial guess vectors are not necessarily orthonormal - - -`davidson_diag_mrcc `_ - Davidson diagonalization. - .br - dets_in : bitmasks corresponding to determinants - .br - u_in : guess coefficients on the various states. Overwritten - on exit - .br - dim_in : leftmost dimension of u_in - .br - sze : Number of determinants - .br - N_st : Number of eigenstates - .br - iunit : Unit number for the I/O - .br - Initial guess vectors are not necessarily orthonormal - - -`delta_ii `_ - Dressing matrix in N_det basis - - -`delta_ij `_ - Dressing matrix in N_det basis - - -`diagonalize_ci_dressed `_ - Replace the coefficients of the CI states by the coefficients of the - eigenstates of the CI matrix - - -`get_excitation_operators_for_one_ref `_ - This subroutine provides all the amplitudes and excitation operators - that one needs to go from the reference to the non reference wave function - you enter with det_ref that is a reference determinant - .br - N_connect_ref is the number of determinants belonging to psi_non_ref - that are connected to det_ref. - .br - amplitudes_phase_less(i) = amplitude phase less t_{I->i} = * lambda_mrcc(i) * phase(I->i) - .br - excitation_operators(:,i) represents the holes and particles that - link the ith connected determinant to det_ref - if :: - excitation_operators(5,i) = 2 :: double excitation alpha - excitation_operators(5,i) = -2 :: double excitation beta - !! excitation_operators(1,i) :: hole 1 - !! excitation_operators(2,i) :: particle 1 - !! excitation_operators(3,i) :: hole 2 - !! excitation_operators(4,i) :: particle 2 - else if :: - excitation_operators(5,i) = 1 :: single excitation alpha - !! excitation_operators(1,i) :: hole 1 - !! excitation_operators(2,i) :: particle 1 - else if :: - excitation_operators(5,i) = -1 :: single excitation beta - !! excitation_operators(3,i) :: hole 1 - !! excitation_operators(4,i) :: particle 1 - else if :: - !! excitation_operators(5,i) = 0 :: double excitation alpha/beta - !! excitation_operators(1,i) :: hole 1 alpha - !! excitation_operators(2,i) :: particle 1 alpha - !! excitation_operators(3,i) :: hole 2 beta - !! excitation_operators(4,i) :: particle 2 beta - - -`h_matrix_dressed `_ - Dressed H with Delta_ij - - -`h_u_0_mrcc `_ - Computes v_0 = H|u_0> - .br - n : number of determinants - .br - H_jj : array of - - -`lambda_mrcc `_ - cm/ or perturbative 1/Delta_E(m) - - -`lambda_pert `_ - cm/ or perturbative 1/Delta_E(m) - - -`mrcc_dress `_ - Undocumented - - -`mrcc_iterations `_ - Undocumented - - -`run_mrcc `_ - Undocumented - - -`set_generators_bitmasks_as_holes_and_particles `_ - Undocumented - diff --git a/plugins/MRCC_Utils_new/davidson.irp.f b/plugins/MRCC_Utils_new/davidson.irp.f deleted file mode 100644 index 0c7bebbd..00000000 --- a/plugins/MRCC_Utils_new/davidson.irp.f +++ /dev/null @@ -1,430 +0,0 @@ -subroutine davidson_diag_mrcc(dets_in,u_in,energies,dim_in,sze,N_st,Nint,iunit,istate) - use bitmasks - implicit none - BEGIN_DOC - ! Davidson diagonalization. - ! - ! 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 - ! - ! iunit : Unit number for the I/O - ! - ! Initial guess vectors are not necessarily orthonormal - END_DOC - integer, intent(in) :: dim_in, sze, N_st, Nint, iunit, istate - 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, allocatable :: H_jj(:) - - double precision :: diag_h_mat_elem - integer :: i - ASSERT (N_st > 0) - ASSERT (sze > 0) - ASSERT (Nint > 0) - ASSERT (Nint == N_int) - PROVIDE mo_bielec_integrals_in_map - allocate(H_jj(sze)) - - !$OMP PARALLEL DEFAULT(NONE) & - !$OMP SHARED(sze,H_jj,N_det_ref,dets_in,Nint,istate,delta_ii,idx_ref) & - !$OMP PRIVATE(i) - !$OMP DO SCHEDULE(guided) - do i=1,sze - H_jj(i) = diag_h_mat_elem(dets_in(1,1,i),Nint) - enddo - !$OMP END DO - !$OMP DO SCHEDULE(guided) - do i=1,N_det_ref - H_jj(idx_ref(i)) += delta_ii(i,istate) - enddo - !$OMP END DO - !$OMP END PARALLEL - - call davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iunit,istate) - deallocate (H_jj) -end - -subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iunit,istate) - use bitmasks - implicit none - BEGIN_DOC - ! Davidson diagonalization with specific diagonal elements of the H matrix - ! - ! H_jj : specific diagonal H matrix elements to diagonalize de Davidson - ! - ! 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 - ! - ! iunit : Unit for the I/O - ! - ! Initial guess vectors are not necessarily orthonormal - END_DOC - integer, intent(in) :: dim_in, sze, N_st, Nint, istate - integer(bit_kind), intent(in) :: dets_in(Nint,2,sze) - double precision, intent(in) :: H_jj(sze) - integer, intent(in) :: iunit - double precision, intent(inout) :: u_in(dim_in,N_st) - double precision, intent(out) :: energies(N_st) - - 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, allocatable :: kl_pairs(:,:) - integer :: k_pairs, kl - - integer :: iter2 - double precision, allocatable :: W(:,:,:), U(:,:,:), R(:,:) - double precision, allocatable :: y(:,:,:,:), h(:,:,:,:), lambda(:) - double precision :: diag_h_mat_elem - double precision :: residual_norm(N_st) - character*(16384) :: write_buffer - double precision :: to_print(2,N_st) - double precision :: cpu, wall - - PROVIDE det_connections - - call write_time(iunit) - call wall_time(wall) - call cpu_time(cpu) - write(iunit,'(A)') '' - write(iunit,'(A)') 'Davidson Diagonalization' - write(iunit,'(A)') '------------------------' - write(iunit,'(A)') '' - call write_int(iunit,N_st,'Number of states') - call write_int(iunit,sze,'Number of determinants') - write(iunit,'(A)') '' - write_buffer = '===== ' - do i=1,N_st - write_buffer = trim(write_buffer)//' ================ ================' - enddo - write(iunit,'(A)') trim(write_buffer) - write_buffer = ' Iter' - do i=1,N_st - write_buffer = trim(write_buffer)//' Energy Residual' - enddo - write(iunit,'(A)') trim(write_buffer) - write_buffer = '===== ' - do i=1,N_st - write_buffer = trim(write_buffer)//' ================ ================' - enddo - write(iunit,'(A)') trim(write_buffer) - - allocate( & - kl_pairs(2,N_st*(N_st+1)/2), & - W(sze,N_st,davidson_sze_max), & - U(sze,N_st,davidson_sze_max), & - R(sze,N_st), & - h(N_st,davidson_sze_max,N_st,davidson_sze_max), & - y(N_st,davidson_sze_max,N_st,davidson_sze_max), & - lambda(N_st*davidson_sze_max)) - - ASSERT (N_st > 0) - ASSERT (sze > 0) - ASSERT (Nint > 0) - ASSERT (Nint == N_int) - - ! Initialization - ! ============== - - k_pairs=0 - do l=1,N_st - do k=1,l - k_pairs+=1 - kl_pairs(1,k_pairs) = k - kl_pairs(2,k_pairs) = l - enddo - enddo - - !$OMP PARALLEL DEFAULT(NONE) & - !$OMP SHARED(U,sze,N_st,overlap,kl_pairs,k_pairs, & - !$OMP Nint,dets_in,u_in) & - !$OMP PRIVATE(k,l,kl,i) - - - ! Orthonormalize initial guess - ! ============================ - - !$OMP DO - do kl=1,k_pairs - k = kl_pairs(1,kl) - l = kl_pairs(2,kl) - if (k/=l) then - overlap(k,l) = u_dot_v(U_in(1,k),U_in(1,l),sze) - overlap(l,k) = overlap(k,l) - else - overlap(k,k) = u_dot_u(U_in(1,k),sze) - endif - enddo - !$OMP END DO - !$OMP END PARALLEL - - call ortho_lowdin(overlap,size(overlap,1),N_st,U_in,size(U_in,1),sze) - - ! Davidson iterations - ! =================== - - converged = .False. - - do while (.not.converged) - - !$OMP PARALLEL DEFAULT(NONE) & - !$OMP PRIVATE(k,i) SHARED(U,u_in,sze,N_st) - do k=1,N_st - !$OMP DO - do i=1,sze - U(i,k,1) = u_in(i,k) - enddo - !$OMP END DO - enddo - !$OMP END PARALLEL - - do iter=1,davidson_sze_max-1 - - ! Compute W_k = H |u_k> - ! ---------------------- - - do k=1,N_st - call H_u_0_mrcc(W(1,k,iter),U(1,k,iter),H_jj,sze,dets_in,Nint,istate) - enddo - - ! Compute h_kl = = - ! ------------------------------------------- - - do l=1,N_st - do k=1,N_st - do iter2=1,iter-1 - h(k,iter2,l,iter) = u_dot_v(U(1,k,iter2),W(1,l,iter),sze) - h(k,iter,l,iter2) = h(k,iter2,l,iter) - enddo - enddo - do k=1,l - h(k,iter,l,iter) = u_dot_v(U(1,k,iter),W(1,l,iter),sze) - h(l,iter,k,iter) = h(k,iter,l,iter) - enddo - enddo - - !DEBUG H MATRIX - !do i=1,iter - ! print '(10(x,F16.10))', h(1,i,1,1:i) - !enddo - !print *, '' - !END - - ! Diagonalize h - ! ------------- - call lapack_diag(lambda,y,h,N_st*davidson_sze_max,N_st*iter) - - ! Express eigenvectors of h in the determinant basis - ! -------------------------------------------------- - - do k=1,N_st - do i=1,sze - U(i,k,iter+1) = 0.d0 - W(i,k,iter+1) = 0.d0 - do l=1,N_st - do iter2=1,iter - U(i,k,iter+1) = U(i,k,iter+1) + U(i,l,iter2)*y(l,iter2,k,1) - W(i,k,iter+1) = W(i,k,iter+1) + W(i,l,iter2)*y(l,iter2,k,1) - enddo - enddo - enddo - enddo - - ! Compute residual vector - ! ----------------------- - - do k=1,N_st - do i=1,sze - R(i,k) = lambda(k) * U(i,k,iter+1) - W(i,k,iter+1) - enddo - residual_norm(k) = u_dot_u(R(1,k),sze) - to_print(1,k) = lambda(k) + nuclear_repulsion - to_print(2,k) = residual_norm(k) - enddo - - write(iunit,'(X,I3,X,100(X,F16.10,X,E16.6))'), iter, to_print(:,1:N_st) - call davidson_converged(lambda,residual_norm,wall,iter,cpu,N_st,converged) - if (converged) then - exit - endif - - - ! Davidson step - ! ------------- - - do k=1,N_st - do i=1,sze - U(i,k,iter+1) = -1.d0/max(H_jj(i) - lambda(k),1.d-2) * R(i,k) - enddo - enddo - - ! Gram-Schmidt - ! ------------ - - double precision :: c - do k=1,N_st - do iter2=1,iter - do l=1,N_st - c = u_dot_v(U(1,k,iter+1),U(1,l,iter2),sze) - do i=1,sze - U(i,k,iter+1) -= c * U(i,l,iter2) - enddo - enddo - enddo - do l=1,k-1 - c = u_dot_v(U(1,k,iter+1),U(1,l,iter+1),sze) - do i=1,sze - U(i,k,iter+1) -= c * U(i,l,iter+1) - enddo - enddo - call normalize( U(1,k,iter+1), sze ) - enddo - - !DEBUG : CHECK OVERLAP - !print *, '===' - !do k=1,iter+1 - ! do l=1,k - ! c = u_dot_v(U(1,1,k),U(1,1,l),sze) - ! print *, k,l, c - ! enddo - !enddo - !print *, '===' - !pause - !END DEBUG - - - enddo - - if (.not.converged) then - iter = davidson_sze_max-1 - endif - - ! Re-contract to u_in - ! ----------- - - do k=1,N_st - energies(k) = lambda(k) - do i=1,sze - u_in(i,k) = 0.d0 - do iter2=1,iter - do l=1,N_st - u_in(i,k) += U(i,l,iter2)*y(l,iter2,k,1) - enddo - enddo - enddo - enddo - - enddo - - write_buffer = '===== ' - do i=1,N_st - write_buffer = trim(write_buffer)//' ================ ================' - enddo - write(iunit,'(A)') trim(write_buffer) - write(iunit,'(A)') '' - call write_time(iunit) - - deallocate ( & - kl_pairs, & - W, & - U, & - R, & - h, & - y, & - lambda & - ) - abort_here = abort_all -end - -subroutine H_u_0_mrcc(v_0,u_0,H_jj,n,keys_tmp,Nint,istate) - use bitmasks - implicit none - BEGIN_DOC - ! Computes v_0 = H|u_0> - ! - ! n : number of determinants - ! - ! H_jj : array of - END_DOC - integer, intent(in) :: n,Nint,istate - double precision, intent(out) :: v_0(n) - double precision, intent(in) :: u_0(n) - double precision, intent(in) :: H_jj(n) - integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n) - integer, allocatable :: idx(:) - double precision :: hij - double precision, allocatable :: vt(:) - integer :: i,j,k,l, jj,ii - integer :: i0, j0 - ASSERT (Nint > 0) - ASSERT (Nint == N_int) - ASSERT (n>0) - PROVIDE ref_bitmask_energy delta_ij - integer, parameter :: block_size = 157 - !$OMP PARALLEL DEFAULT(NONE) & - !$OMP PRIVATE(i,hij,j,k,idx,jj,ii,vt) & - !$OMP SHARED(n_det_ref,n_det_non_ref,idx_ref,idx_non_ref,n,H_jj,u_0,keys_tmp,Nint,v_0,istate,delta_ij) - !$OMP DO SCHEDULE(static) - do i=1,n - v_0(i) = H_jj(i) * u_0(i) - enddo - !$OMP END DO - allocate(idx(0:n), vt(n)) - Vt = 0.d0 - !$OMP DO SCHEDULE(guided) - do i=1,n - idx(0) = i - call filter_connected_davidson(keys_tmp,keys_tmp(1,1,i),Nint,i-1,idx) - do jj=1,idx(0) - j = idx(jj) - if ( (dabs(u_0(j)) > 1.d-7).or.((dabs(u_0(i)) > 1.d-7)) ) then - call i_H_j(keys_tmp(1,1,j),keys_tmp(1,1,i),Nint,hij) - hij = hij - vt (i) = vt (i) + hij*u_0(j) - vt (j) = vt (j) + hij*u_0(i) - endif - enddo - enddo - !$OMP END DO - - !$OMP DO SCHEDULE(guided) - do ii=1,n_det_ref - i = idx_ref(ii) - do jj = 1, n_det_non_ref - j = idx_non_ref(jj) - vt (i) = vt (i) + delta_ij(ii,jj,istate)*u_0(j) - vt (j) = vt (j) + delta_ij(ii,jj,istate)*u_0(i) - enddo - enddo - !$OMP END DO - !$OMP CRITICAL - do i=1,n - v_0(i) = v_0(i) + vt(i) - enddo - !$OMP END CRITICAL - deallocate(idx,vt) - !$OMP END PARALLEL -end - - diff --git a/plugins/MRCC_Utils_new/mrcc_amplitudes.irp.f b/plugins/MRCC_Utils_new/mrcc_amplitudes.irp.f deleted file mode 100644 index 6746bee1..00000000 --- a/plugins/MRCC_Utils_new/mrcc_amplitudes.irp.f +++ /dev/null @@ -1,93 +0,0 @@ -subroutine get_excitation_operators_for_one_ref(det_ref,i_state,ndetnonref,N_connect_ref,excitation_operators,amplitudes_phase_less,index_connected) - use bitmasks - implicit none - integer(bit_kind), intent(in) :: det_ref(N_int,2) - integer, intent(in) :: i_state,ndetnonref - integer*2, intent(out) :: excitation_operators(5,ndetnonref) - integer, intent(out) :: index_connected(ndetnonref) - integer, intent(out) :: N_connect_ref - double precision, intent(out) :: amplitudes_phase_less(ndetnonref) - - integer :: i,j,k,l,degree,h1,p1,h2,p2,s1,s2 - integer :: exc(0:2,2,2) - double precision :: phase,hij - BEGIN_DOC - ! This subroutine provides all the amplitudes and excitation operators - ! that one needs to go from the reference to the non reference wave function - ! you enter with det_ref that is a reference determinant - ! - ! N_connect_ref is the number of determinants belonging to psi_non_ref - ! that are connected to det_ref. - ! - ! amplitudes_phase_less(i) = amplitude phase less t_{I->i} = * lambda_mrcc(i) * phase(I->i) - ! - ! excitation_operators(:,i) represents the holes and particles that - ! link the ith connected determinant to det_ref - ! if :: - ! excitation_operators(5,i) = 2 :: double excitation alpha - ! excitation_operators(5,i) = -2 :: double excitation beta - !!! excitation_operators(1,i) :: hole 1 - !!! excitation_operators(2,i) :: particle 1 - !!! excitation_operators(3,i) :: hole 2 - !!! excitation_operators(4,i) :: particle 2 - ! else if :: - ! excitation_operators(5,i) = 1 :: single excitation alpha - !!! excitation_operators(1,i) :: hole 1 - !!! excitation_operators(2,i) :: particle 1 - ! else if :: - ! excitation_operators(5,i) = -1 :: single excitation beta - !!! excitation_operators(3,i) :: hole 1 - !!! excitation_operators(4,i) :: particle 1 - ! else if :: - !!! excitation_operators(5,i) = 0 :: double excitation alpha/beta - !!! excitation_operators(1,i) :: hole 1 alpha - !!! excitation_operators(2,i) :: particle 1 alpha - !!! excitation_operators(3,i) :: hole 2 beta - !!! excitation_operators(4,i) :: particle 2 beta - END_DOC - N_connect_ref = 0 - do i = 1, ndetnonref - call i_H_j_phase_out(det_ref,psi_non_ref(1,1,i),N_int,hij,phase,exc,degree) - if (dabs(hij) <= mo_integrals_threshold) then - cycle - endif - N_connect_ref +=1 - index_connected(N_connect_ref) = i - call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) - amplitudes_phase_less(N_connect_ref) = hij * lambda_mrcc(i_state,i) !*phase - - if (degree==2) then - - excitation_operators(1,N_connect_ref) = h1 - excitation_operators(2,N_connect_ref) = p1 - excitation_operators(3,N_connect_ref) = h2 - excitation_operators(4,N_connect_ref) = p2 - if(s1==s2.and.s1==1)then ! double alpha - excitation_operators(5,N_connect_ref) = 2 - elseif(s1==s2.and.s1==2)then ! double beta - excitation_operators(5,N_connect_ref) = -2 - else ! double alpha/beta - excitation_operators(5,N_connect_ref) = 0 - endif - - else if(degree==1) then - - if(s1==1)then ! mono alpha - excitation_operators(5,N_connect_ref) = 1 - excitation_operators(1,N_connect_ref) = h1 - excitation_operators(2,N_connect_ref) = p1 - else ! mono beta - excitation_operators(5,N_connect_ref) = -1 - excitation_operators(3,N_connect_ref) = h1 - excitation_operators(4,N_connect_ref) = p1 - endif - - else - - N_connect_ref-=1 - - endif - - enddo - -end diff --git a/plugins/MRCC_Utils_new/mrcc_dress.irp.f b/plugins/MRCC_Utils_new/mrcc_dress.irp.f deleted file mode 100644 index ee998995..00000000 --- a/plugins/MRCC_Utils_new/mrcc_dress.irp.f +++ /dev/null @@ -1,183 +0,0 @@ -subroutine mrcc_dress(ndetref,ndetnonref,nstates,delta_ij_,delta_ii_) - use bitmasks - implicit none - integer, intent(in) :: ndetref,nstates,ndetnonref - double precision, intent(inout) :: delta_ii_(ndetref,nstates),delta_ij_(ndetref,ndetnonref,nstates) - integer :: i,j,k,l,m - integer :: i_state - integer :: N_connect_ref - integer*2,allocatable :: excitation_operators(:,:) - double precision, allocatable :: amplitudes_phase_less(:) - double precision, allocatable :: coef_test(:) - integer(bit_kind), allocatable :: key_test(:,:) - integer, allocatable :: index_connected(:) - integer :: i_hole,i_particle,ispin,i_ok,connected_to_ref,index_wf - integer, allocatable :: idx_vector(:) - double precision :: phase_ij - double precision :: dij,phase_la - double precision :: hij,phase - integer :: exc(0:2,2,2),degree - logical :: is_in_wavefunction - double precision, allocatable :: delta_ij_tmp(:,:,:), delta_ii_tmp(:,:) - logical, external :: is_in_psi_ref - - i_state = 1 - allocate(excitation_operators(5,N_det_non_ref)) - allocate(amplitudes_phase_less(N_det_non_ref)) - allocate(index_connected(N_det_non_ref)) - - !$OMP PARALLEL DEFAULT(NONE) & - !$OMP SHARED(N_det_ref, N_det_non_ref, psi_ref, i_state, & - !$OMP N_connect_ref,index_connected,psi_non_ref, & - !$OMP excitation_operators,amplitudes_phase_less, & - !$OMP psi_non_ref_coef,N_int,lambda_mrcc, & - !$OMP delta_ii_,delta_ij_,psi_ref_coef,nstates, & - !$OMP mo_integrals_threshold,idx_non_ref_rev) & - !$OMP PRIVATE(i,j,k,l,hil,phase_il,exc,degree,t_il, & - !$OMP key_test,i_ok,phase_la,hij,phase_ij,m, & - !$OMP dij,idx_vector,delta_ij_tmp, & - !$OMP delta_ii_tmp,phase) - allocate(idx_vector(0:N_det_non_ref)) - allocate(key_test(N_int,2)) - allocate(delta_ij_tmp(size(delta_ij_,1),size(delta_ij_,2),nstates)) - allocate(delta_ii_tmp(size(delta_ij_,1),nstates)) - delta_ij_tmp = 0.d0 - delta_ii_tmp = 0.d0 - - do i = 1, N_det_ref - !$OMP SINGLE - call get_excitation_operators_for_one_ref(psi_ref(1,1,i),i_state,N_det_non_ref,N_connect_ref,excitation_operators,amplitudes_phase_less,index_connected) - print*,'N_connect_ref =',N_connect_ref - print*,'N_det_non_ref =',N_det_non_ref - !$OMP END SINGLE - !$OMP BARRIER - - !$OMP DO SCHEDULE(dynamic) - do l = 1, N_det_non_ref -! print *, l, '/', N_det_non_ref - double precision :: t_il,phase_il,hil - call i_H_j_phase_out(psi_ref(1,1,i),psi_non_ref(1,1,l),N_int,hil,phase_il,exc,degree) - t_il = hil * lambda_mrcc(i_state,l) - if (dabs(t_il) < mo_integrals_threshold) then - cycle - endif - ! loop on the non ref determinants - - do j = 1, N_connect_ref - ! loop on the excitation operators linked to i - - do k = 1, N_int - key_test(k,1) = psi_non_ref(k,1,l) - key_test(k,2) = psi_non_ref(k,2,l) - enddo - - ! we apply the excitation operator T_I->j - call apply_excitation_operator(key_test,excitation_operators(1,j),i_ok) - if(i_ok.ne.1)cycle - - ! we check if such determinant is already in the wave function - if(is_in_wavefunction(key_test,N_int))cycle - - ! we get the phase for psi_non_ref(l) -> T_I->j |psi_non_ref(l)> - call get_excitation(psi_non_ref(1,1,l),key_test,exc,degree,phase_la,N_int) - - ! we get the phase T_I->j - call i_H_j_phase_out(psi_ref(1,1,i),psi_non_ref(1,1,index_connected(j)),N_int,hij,phase_ij,exc,degree) - - ! we compute the contribution to the coef of key_test - dij = t_il * hij * phase_la *phase_ij *lambda_mrcc(i_state,index_connected(j)) * 0.5d0 - if (dabs(dij) < mo_integrals_threshold) then - cycle - endif - - ! we compute the interaction of such determinant with all the non_ref dets - call filter_connected(psi_non_ref,key_test,N_int,N_det_non_ref,idx_vector) - - do k = 1, idx_vector(0) - m = idx_vector(k) - call i_H_j_phase_out(key_test,psi_non_ref(1,1,m),N_int,hij,phase,exc,degree) - delta_ij_tmp(i,m,i_state) += hij * dij - enddo - - - enddo - - if(dabs(psi_ref_coef(i,i_state)).le.5.d-5) then - delta_ii_tmp(i,i_state) -= & - delta_ij_tmp(i,l,i_state) * psi_non_ref_coef(l,i_state) & - / psi_ref_coef(i,i_state) - endif - - enddo - !$OMP END DO - enddo - - !$OMP CRITICAL - delta_ij_ = delta_ij_ + delta_ij_tmp - delta_ii_ = delta_ii_ + delta_ii_tmp - !$OMP END CRITICAL - - deallocate(delta_ii_tmp,delta_ij_tmp) - deallocate(idx_vector) - deallocate(key_test) - !$OMP END PARALLEL - - deallocate(excitation_operators) - deallocate(amplitudes_phase_less) - -end - - - -subroutine apply_excitation_operator(key_in,excitation_operator,i_ok) - use bitmasks - implicit none - integer(bit_kind), intent(inout) :: key_in - integer, intent (out) :: i_ok - integer*2 :: excitation_operator(5) - integer :: i_particle,i_hole,ispin - ! Do excitation - if(excitation_operator(5)==1)then ! mono alpha - i_hole = excitation_operator(1) - i_particle = excitation_operator(2) - ispin = 1 - call do_mono_excitation(key_in,i_hole,i_particle,ispin,i_ok) - else if (excitation_operator(5)==-1)then ! mono beta - i_hole = excitation_operator(3) - i_particle = excitation_operator(4) - ispin = 2 - call do_mono_excitation(key_in,i_hole,i_particle,ispin,i_ok) - else if (excitation_operator(5) == -2 )then ! double beta - i_hole = excitation_operator(1) - i_particle = excitation_operator(2) - ispin = 2 - call do_mono_excitation(key_in,i_hole,i_particle,ispin,i_ok) - if(i_ok.ne.1)return - i_hole = excitation_operator(3) - i_particle = excitation_operator(4) - ispin = 2 - call do_mono_excitation(key_in,i_hole,i_particle,ispin,i_ok) - - else if (excitation_operator(5) == 2 )then ! double alpha - i_hole = excitation_operator(1) - i_particle = excitation_operator(2) - ispin = 1 - call do_mono_excitation(key_in,i_hole,i_particle,ispin,i_ok) - if(i_ok.ne.1)return - i_hole = excitation_operator(3) - i_particle = excitation_operator(4) - ispin = 1 - call do_mono_excitation(key_in,i_hole,i_particle,ispin,i_ok) - - else if (excitation_operator(5) == 0 )then ! double alpha/alpha - i_hole = excitation_operator(1) - i_particle = excitation_operator(2) - ispin = 1 - call do_mono_excitation(key_in,i_hole,i_particle,ispin,i_ok) - if(i_ok.ne.1)return - i_hole = excitation_operator(3) - i_particle = excitation_operator(4) - ispin = 2 - call do_mono_excitation(key_in,i_hole,i_particle,ispin,i_ok) - endif -end diff --git a/plugins/MRCC_Utils_new/mrcc_general.irp.f b/plugins/MRCC_Utils_new/mrcc_general.irp.f deleted file mode 100644 index 245fcb05..00000000 --- a/plugins/MRCC_Utils_new/mrcc_general.irp.f +++ /dev/null @@ -1,67 +0,0 @@ -subroutine run_mrcc - implicit none - call set_generators_bitmasks_as_holes_and_particles - call mrcc_iterations -end - -subroutine mrcc_iterations - implicit none - - integer :: i,j - - double precision :: E_new, E_old, delta_e - integer :: iteration - E_new = 0.d0 - delta_E = 1.d0 - iteration = 0 - do while (delta_E > 1.d-8) - iteration += 1 - print *, '===========================' - print *, 'MRCC Iteration', iteration - print *, '===========================' - print *, '' - E_old = sum(ci_energy_dressed) - call write_double(6,ci_energy_dressed(1),"MRCC energy") - call diagonalize_ci_dressed - E_new = sum(ci_energy_dressed) - delta_E = dabs(E_new - E_old) -! stop - if (iteration > 200) then - exit - endif - enddo - call write_double(6,ci_energy_dressed(1),"Final MRCC energy") - call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1)) - call save_wavefunction - -end - -subroutine set_generators_bitmasks_as_holes_and_particles - implicit none - integer :: i,k - do k = 1, N_generators_bitmask - do i = 1, N_int - ! Pure single part - generators_bitmask(i,1,1,k) = holes_operators(i,1) ! holes for pure single exc alpha - generators_bitmask(i,1,2,k) = particles_operators(i,1) ! particles for pure single exc alpha - generators_bitmask(i,2,1,k) = holes_operators(i,2) ! holes for pure single exc beta - generators_bitmask(i,2,2,k) = particles_operators(i,2) ! particles for pure single exc beta - - ! Double excitation - generators_bitmask(i,1,3,k) = holes_operators(i,1) ! holes for first single exc alpha - generators_bitmask(i,1,4,k) = particles_operators(i,1) ! particles for first single exc alpha - generators_bitmask(i,2,3,k) = holes_operators(i,2) ! holes for first single exc beta - generators_bitmask(i,2,4,k) = particles_operators(i,2) ! particles for first single exc beta - - generators_bitmask(i,1,5,k) = holes_operators(i,1) ! holes for second single exc alpha - generators_bitmask(i,1,6,k) = particles_operators(i,1) ! particles for second single exc alpha - generators_bitmask(i,2,5,k) = holes_operators(i,2) ! holes for second single exc beta - generators_bitmask(i,2,6,k) = particles_operators(i,2) ! particles for second single exc beta - - enddo - enddo - touch generators_bitmask - - - -end diff --git a/plugins/MRCC_Utils_new/mrcc_utils.irp.f b/plugins/MRCC_Utils_new/mrcc_utils.irp.f deleted file mode 100644 index d97696e5..00000000 --- a/plugins/MRCC_Utils_new/mrcc_utils.irp.f +++ /dev/null @@ -1,179 +0,0 @@ - BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ] -&BEGIN_PROVIDER [ double precision, lambda_pert, (N_states,psi_det_size) ] - implicit none - BEGIN_DOC - ! cm/ or perturbative 1/Delta_E(m) - END_DOC - integer :: i,k - double precision :: ihpsi(N_states), hii - integer :: i_ok - i_ok = 0 - - do i=1,N_det_non_ref - call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,& - size(psi_ref_coef,1), n_states, ihpsi) - call i_h_j(psi_non_ref(1,1,i),psi_non_ref(1,1,i),N_int,hii) - do k=1,N_states - lambda_pert(k,i) = 1.d0 / (psi_ref_energy_diagonalized(k)-hii) - if (dabs(ihpsi(k)).le.1.d-3) then - i_ok +=1 - lambda_mrcc(k,i) = lambda_pert(k,i) - else - lambda_mrcc(k,i) = psi_non_ref_coef(i,k)/ihpsi(k) - endif - enddo - enddo - print*,'N_det_non_ref = ',N_det_non_ref - print*,'Number of Perturbatively treated determinants = ',i_ok - print*,'psi_coef_ref_ratio = ',psi_ref_coef(2,1)/psi_ref_coef(1,1) - -END_PROVIDER - - - - -!BEGIN_PROVIDER [ double precision, delta_ij_non_ref, (N_det_non_ref, N_det_non_ref,N_states) ] -!implicit none -!BEGIN_DOC -!! Dressing matrix in SD basis -!END_DOC -!delta_ij_non_ref = 0.d0 -!call H_apply_mrcc_simple(delta_ij_non_ref,N_det_non_ref) -!END_PROVIDER - - BEGIN_PROVIDER [ double precision, delta_ij, (N_det_ref,N_det_non_ref,N_states) ] -&BEGIN_PROVIDER [ double precision, delta_ii, (N_det_ref,N_states) ] - implicit none - BEGIN_DOC - ! Dressing matrix in N_det basis - END_DOC - integer :: i,j,m - delta_ij = 0.d0 - delta_ii = 0.d0 - call mrcc_dress(N_det_ref,N_det_non_ref,N_states,delta_ij,delta_ii) - write(33,*)delta_ij - write(34,*)delta_ii -END_PROVIDER - -BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det,N_states) ] - implicit none - BEGIN_DOC - ! Dressed H with Delta_ij - END_DOC - integer :: i, j,istate,ii,jj - do istate = 1,N_states - do j=1,N_det - do i=1,N_det - h_matrix_dressed(i,j,istate) = h_matrix_all_dets(i,j) - enddo - enddo - do ii = 1, N_det_ref - i =idx_ref(ii) - h_matrix_dressed(i,i,istate) += delta_ii(ii,istate) - do jj = 1, N_det_non_ref - j =idx_non_ref(jj) - h_matrix_dressed(i,j,istate) += delta_ij(ii,jj,istate) - h_matrix_dressed(j,i,istate) += delta_ij(ii,jj,istate) - enddo - enddo - enddo -END_PROVIDER - - - BEGIN_PROVIDER [ double precision, CI_electronic_energy_dressed, (N_states_diag) ] -&BEGIN_PROVIDER [ double precision, CI_eigenvectors_dressed, (N_det,N_states_diag) ] -&BEGIN_PROVIDER [ double precision, CI_eigenvectors_s2_dressed, (N_states_diag) ] - implicit none - BEGIN_DOC - ! Eigenvectors/values of the CI matrix - END_DOC - integer :: i,j - - do j=1,N_states_diag - do i=1,N_det - CI_eigenvectors_dressed(i,j) = psi_coef(i,j) - enddo - enddo - - if (diag_algorithm == "Davidson") then - - integer :: istate - istate = 1 - call davidson_diag_mrcc(psi_det,CI_eigenvectors_dressed,CI_electronic_energy_dressed,& - size(CI_eigenvectors_dressed,1),N_det,N_states_diag,N_int,output_determinants,istate) - - else if (diag_algorithm == "Lapack") then - - double precision, allocatable :: eigenvectors(:,:), eigenvalues(:) - allocate (eigenvectors(size(H_matrix_dressed,1),N_det)) - allocate (eigenvalues(N_det)) - call lapack_diag(eigenvalues,eigenvectors, & - H_matrix_dressed,size(H_matrix_dressed,1),N_det) - CI_electronic_energy_dressed(:) = 0.d0 - do i=1,N_det - CI_eigenvectors_dressed(i,1) = eigenvectors(i,1) - enddo - integer :: i_state - double precision :: s2 - i_state = 0 - if (s2_eig) then - do j=1,N_det - call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2) - if(dabs(s2-expected_s2).le.0.3d0)then - i_state += 1 - do i=1,N_det - CI_eigenvectors_dressed(i,i_state) = eigenvectors(i,j) - enddo - CI_electronic_energy_dressed(i_state) = eigenvalues(j) - CI_eigenvectors_s2_dressed(i_state) = s2 - endif - if (i_state.ge.N_states_diag) then - exit - endif - enddo - else - do j=1,N_states_diag - call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2) - i_state += 1 - do i=1,N_det - CI_eigenvectors_dressed(i,i_state) = eigenvectors(i,j) - enddo - CI_electronic_energy_dressed(i_state) = eigenvalues(j) - CI_eigenvectors_s2_dressed(i_state) = s2 - enddo - endif - deallocate(eigenvectors,eigenvalues) - endif - -END_PROVIDER - -BEGIN_PROVIDER [ double precision, CI_energy_dressed, (N_states_diag) ] - implicit none - BEGIN_DOC - ! N_states lowest eigenvalues of the dressed CI matrix - END_DOC - - integer :: j - character*(8) :: st - call write_time(output_determinants) - do j=1,N_states_diag - CI_energy_dressed(j) = CI_electronic_energy_dressed(j) + nuclear_repulsion - enddo - -END_PROVIDER - -subroutine diagonalize_CI_dressed - implicit none - BEGIN_DOC -! Replace the coefficients of the CI states by the coefficients of the -! eigenstates of the CI matrix - END_DOC - integer :: i,j - do j=1,N_states_diag - do i=1,N_det - psi_coef(i,j) = CI_eigenvectors_dressed(i,j) - enddo - enddo - SOFT_TOUCH psi_coef - -end diff --git a/plugins/MRCC_Utils_new/tree_dependency.png b/plugins/MRCC_Utils_new/tree_dependency.png deleted file mode 100644 index 500e5d43..00000000 Binary files a/plugins/MRCC_Utils_new/tree_dependency.png and /dev/null differ diff --git a/plugins/QmcChem/dressed_dmc.irp.f b/plugins/QmcChem/dressed_dmc.irp.f index 0a48e871..1437ed31 100644 --- a/plugins/QmcChem/dressed_dmc.irp.f +++ b/plugins/QmcChem/dressed_dmc.irp.f @@ -57,7 +57,7 @@ program dressed_dmc enddo - call davidson_diag_hjj(psi_det,psi_coef,H_jj,energies,size(psi_coef,1),N_det,N_states,N_int,6) + call davidson_diag_hjj(psi_det,psi_coef,H_jj,energies,size(psi_coef,1),N_det,N_states,N_states_diag,,N_int,6) call save_wavefunction call write_spindeterminants diff --git a/plugins/mrcepa0/EZFIO.cfg b/plugins/mrcepa0/EZFIO.cfg index 903e6e95..7580f028 100644 --- a/plugins/mrcepa0/EZFIO.cfg +++ b/plugins/mrcepa0/EZFIO.cfg @@ -1,5 +1,5 @@ [lambda_type] type: Positive_int -doc: lambda type ( 0 = none, 1 = last version ) +doc: lambda type interface: ezfio,provider,ocaml default: 0 diff --git a/scripts/compilation/qp_create_ninja.py b/scripts/compilation/qp_create_ninja.py index e2acbf7f..b5ae96e5 100755 --- a/scripts/compilation/qp_create_ninja.py +++ b/scripts/compilation/qp_create_ninja.py @@ -796,7 +796,7 @@ def create_build_ninja_global(): l_string += ["build dummy_target: update_build_ninja_root", - "build ocaml_target: make_ocaml", + "build ocaml_target: make_ocaml | dummy_target", "", "build all: make_all dummy_target ocaml_target", "default all", diff --git a/src/Davidson/diagonalization.irp.f b/src/Davidson/diagonalization.irp.f index 8ece1c5c..21b13113 100644 --- a/src/Davidson/diagonalization.irp.f +++ b/src/Davidson/diagonalization.irp.f @@ -1,4 +1,4 @@ -subroutine davidson_diag(dets_in,u_in,energies,dim_in,sze,N_st,Nint,iunit) +subroutine davidson_diag(dets_in,u_in,energies,dim_in,sze,N_st,N_st_diag,Nint,iunit) use bitmasks implicit none BEGIN_DOC @@ -19,9 +19,9 @@ subroutine davidson_diag(dets_in,u_in,energies,dim_in,sze,N_st,Nint,iunit) ! ! Initial guess vectors are not necessarily orthonormal END_DOC - integer, intent(in) :: dim_in, sze, N_st, Nint, iunit + integer, intent(in) :: dim_in, sze, N_st, N_st_diag, Nint, iunit integer(bit_kind), intent(in) :: dets_in(Nint,2,sze) - double precision, intent(inout) :: u_in(dim_in,N_st) + double precision, intent(inout) :: u_in(dim_in,N_st_diag) double precision, intent(out) :: energies(N_st) double precision, allocatable :: H_jj(:) @@ -44,7 +44,7 @@ subroutine davidson_diag(dets_in,u_in,energies,dim_in,sze,N_st,Nint,iunit) !$OMP END DO !$OMP END PARALLEL - call davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iunit) + call davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag,Nint,iunit) deallocate (H_jj) end @@ -270,7 +270,7 @@ subroutine sort_dets_ab(key, idx, shortcut, N_key, Nint) end subroutine -subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iunit) +subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag,Nint,iunit) use bitmasks implicit none BEGIN_DOC @@ -288,24 +288,26 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun ! sze : Number of determinants ! ! N_st : Number of eigenstates + ! + ! N_st_diag : Number of states in which H is diagonalized ! ! iunit : Unit for the I/O ! ! Initial guess vectors are not necessarily orthonormal END_DOC - integer, intent(in) :: dim_in, sze, N_st, Nint + integer, intent(in) :: dim_in, sze, N_st, N_st_diag, Nint integer(bit_kind), intent(in) :: dets_in(Nint,2,sze) double precision, intent(in) :: H_jj(sze) integer, intent(in) :: iunit - double precision, intent(inout) :: u_in(dim_in,N_st) - double precision, intent(out) :: energies(N_st) + double precision, intent(inout) :: u_in(dim_in,N_st_diag) + double precision, intent(out) :: energies(N_st_diag) integer :: sze_8 integer :: iter integer :: i,j,k,l,m logical :: converged - double precision :: overlap(N_st,N_st) + double precision, allocatable :: overlap(:,:) double precision :: u_dot_v, u_dot_u integer, allocatable :: kl_pairs(:,:) @@ -315,13 +317,14 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun double precision, allocatable :: W(:,:,:), U(:,:,:), R(:,:) double precision, allocatable :: y(:,:,:,:), h(:,:,:,:), lambda(:) double precision :: diag_h_mat_elem - double precision :: residual_norm(N_st) + double precision, allocatable :: residual_norm(:) character*(16384) :: write_buffer double precision :: to_print(2,N_st) double precision :: cpu, wall !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: U, W, R, y, h, lambda + PROVIDE nuclear_repulsion call write_time(iunit) call wall_time(wall) @@ -331,6 +334,7 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun write(iunit,'(A)') '------------------------' write(iunit,'(A)') '' call write_int(iunit,N_st,'Number of states') + call write_int(iunit,N_st_diag,'Number of states in diagonalization') call write_int(iunit,sze,'Number of determinants') write(iunit,'(A)') '' write_buffer = '===== ' @@ -353,70 +357,22 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun sze_8 = align_double(sze) allocate( & - kl_pairs(2,N_st*(N_st+1)/2), & - W(sze_8,N_st,davidson_sze_max), & - U(sze_8,N_st,davidson_sze_max), & - R(sze_8,N_st), & - h(N_st,davidson_sze_max,N_st,davidson_sze_max), & - y(N_st,davidson_sze_max,N_st,davidson_sze_max), & - lambda(N_st*davidson_sze_max)) + kl_pairs(2,N_st_diag*(N_st_diag+1)/2), & + W(sze_8,N_st_diag,davidson_sze_max), & + U(sze_8,N_st_diag,davidson_sze_max), & + R(sze_8,N_st_diag), & + h(N_st_diag,davidson_sze_max,N_st_diag,davidson_sze_max), & + y(N_st_diag,davidson_sze_max,N_st_diag,davidson_sze_max), & + residual_norm(N_st_diag), & + overlap(N_st_diag,N_st_diag), & + lambda(N_st_diag*davidson_sze_max)) ASSERT (N_st > 0) + ASSERT (N_st_diag >= N_st) ASSERT (sze > 0) ASSERT (Nint > 0) ASSERT (Nint == N_int) - ! Initialization - ! ============== - - - if (N_st > 1) then - - k_pairs=0 - do l=1,N_st - do k=1,l - k_pairs+=1 - kl_pairs(1,k_pairs) = k - kl_pairs(2,k_pairs) = l - enddo - enddo - - !$OMP PARALLEL DEFAULT(NONE) & - !$OMP SHARED(U,sze,N_st,overlap,kl_pairs,k_pairs, & - !$OMP Nint,dets_in,u_in) & - !$OMP PRIVATE(k,l,kl) - - - ! Orthonormalize initial guess - ! ============================ - - !$OMP DO - do kl=1,k_pairs - k = kl_pairs(1,kl) - l = kl_pairs(2,kl) - if (k/=l) then - overlap(k,l) = u_dot_v(U_in(1,k),U_in(1,l),sze) - overlap(l,k) = overlap(k,l) - else - overlap(k,k) = u_dot_u(U_in(1,k),sze) - endif - enddo - !$OMP END DO - !$OMP END PARALLEL - - call ortho_lowdin(overlap,size(overlap,1),N_st,U_in,size(U_in,1),sze) - - else - - overlap(1,1) = u_dot_u(U_in(1,1),sze) - double precision :: f - f = 1.d0 / dsqrt(overlap(1,1)) - do i=1,sze - U_in(i,1) = U_in(i,1) * f - enddo - - endif - ! Davidson iterations ! =================== @@ -424,40 +380,57 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun do while (.not.converged) - !$OMP PARALLEL DEFAULT(NONE) & - !$OMP PRIVATE(k,i) SHARED(U,u_in,sze,N_st) do k=1,N_st - !$OMP DO do i=1,sze U(i,k,1) = u_in(i,k) enddo - !$OMP END DO enddo - !$OMP END PARALLEL + + do k=N_st+1,N_st_diag + do i=1,sze + call RANDOM_NUMBER(U(i,k,1)) + U(i,k,1) = U(i,k,1) - 0.5d0 + enddo + enddo + ! Compute overlap matrix + call DGEMM('T','N', N_st_diag, N_st_diag, sze, & + 1.d0, U(1,1,1), size(U,1), U(1,1,1), size(U,1), & + 0.d0, overlap, size(overlap,1)) + + call ortho_lowdin(overlap,size(overlap,1),N_st_diag,U,size(U,1),sze) + do iter=1,davidson_sze_max-1 ! Compute W_k = H |u_k> ! ---------------------- - call H_u_0_nstates(W(1,1,iter),U(1,1,iter),H_jj,sze,dets_in,Nint,N_st,sze_8) + call H_u_0_nstates(W(1,1,iter),U(1,1,iter),H_jj,sze,dets_in,Nint,N_st_diag,sze_8) ! Compute h_kl = = ! ------------------------------------------- - do l=1,N_st - do k=1,N_st + !$OMP PARALLEL PRIVATE(k,l,iter2) SHARED(h,U,W,sze,iter) + do l=1,N_st_diag + !$OMP DO + do k=1,N_st_diag do iter2=1,iter-1 h(k,iter2,l,iter) = u_dot_v(U(1,k,iter2),W(1,l,iter),sze) h(k,iter,l,iter2) = h(k,iter2,l,iter) enddo enddo + !$OMP END DO + enddo + do l=1,N_st_diag + !$OMP DO do k=1,l h(k,iter,l,iter) = u_dot_v(U(1,k,iter),W(1,l,iter),sze) h(l,iter,k,iter) = h(k,iter,l,iter) enddo + !$OMP END DO enddo + !$OMP END PARALLEL !DEBUG H MATRIX !do i=1,iter @@ -468,16 +441,16 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun ! Diagonalize h ! ------------- - call lapack_diag(lambda,y,h,N_st*davidson_sze_max,N_st*iter) + call lapack_diag(lambda,y,h,N_st_diag*davidson_sze_max,N_st_diag*iter) ! Express eigenvectors of h in the determinant basis ! -------------------------------------------------- !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(k,i,l,iter2) & - !$OMP SHARED(U,W,R,y,iter,lambda,N_st,sze,to_print, & - !$OMP residual_norm,nuclear_repulsion) - do k=1,N_st + !$OMP SHARED(U,W,R,y,iter,lambda,N_st_diag,sze,to_print, & + !$OMP residual_norm,nuclear_repulsion,N_st) + do k=1,N_st_diag !$OMP DO do i=1,sze U(i,k,iter+1) = 0.d0 @@ -485,7 +458,7 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun enddo !$OMP END DO do iter2=1,iter - do l=1,N_st + do l=1,N_st_diag !$OMP DO do i=1,sze U(i,k,iter+1) = U(i,k,iter+1) + U(i,l,iter2)*y(l,iter2,k,1) @@ -504,15 +477,17 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun enddo !$OMP END DO !$OMP SINGLE - residual_norm(k) = u_dot_u(R(1,k),sze) - to_print(1,k) = lambda(k) + nuclear_repulsion - to_print(2,k) = residual_norm(k) + if (k <= N_st) then + residual_norm(k) = u_dot_u(R(1,k),sze) + to_print(1,k) = lambda(k) + nuclear_repulsion + to_print(2,k) = residual_norm(k) + endif !$OMP END SINGLE enddo !$OMP END PARALLEL write(iunit,'(X,I3,X,100(X,F16.10,X,E16.6))') iter, to_print(:,1:N_st) - call davidson_converged(lambda,residual_norm,wall,iter,cpu,N_states_diag,converged) + call davidson_converged(lambda,residual_norm,wall,iter,cpu,N_st,converged) if (converged) then exit endif @@ -521,7 +496,7 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun ! Davidson step ! ------------- - do k=1,N_st + do k=1,N_st_diag do i=1,sze U(i,k,iter+1) = -1.d0/max(H_jj(i) - lambda(k),1.d-2) * R(i,k) enddo @@ -531,9 +506,9 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun ! ------------ double precision :: c - do k=1,N_st + do k=1,N_st_diag do iter2=1,iter - do l=1,N_st + do l=1,N_st_diag c = u_dot_v(U(1,k,iter+1),U(1,l,iter2),sze) do i=1,sze U(i,k,iter+1) = U(i,k,iter+1) - c * U(i,l,iter2) @@ -571,12 +546,12 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun ! Re-contract to u_in ! ----------- - do k=1,N_st + do k=1,N_st_diag energies(k) = lambda(k) do i=1,sze u_in(i,k) = 0.d0 do iter2=1,iter - do l=1,N_st + do l=1,N_st_diag u_in(i,k) += U(i,l,iter2)*y(l,iter2,k,1) enddo enddo @@ -595,8 +570,8 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun deallocate ( & kl_pairs, & - W, & - U, & + W, residual_norm, & + U, overlap, & R, & h, & y, & diff --git a/src/Davidson/diagonalize_CI.irp.f b/src/Davidson/diagonalize_CI.irp.f index 291fd0bc..dcf8101d 100644 --- a/src/Davidson/diagonalize_CI.irp.f +++ b/src/Davidson/diagonalize_CI.irp.f @@ -54,7 +54,7 @@ END_PROVIDER if (diag_algorithm == "Davidson") then call davidson_diag(psi_det,CI_eigenvectors,CI_electronic_energy,& - size(CI_eigenvectors,1),N_det,N_states_diag,N_int,output_determinants) + size(CI_eigenvectors,1),N_det,N_states,N_states_diag,N_int,output_determinants) call u_0_S2_u_0(CI_eigenvectors_s2,CI_eigenvectors,N_det,psi_det,N_int,& N_states_diag,size(CI_eigenvectors,1)) diff --git a/src/Davidson/diagonalize_CI_mono.irp.f b/src/Davidson/diagonalize_CI_mono.irp.f index 3f9b94ec..b45d03d6 100644 --- a/src/Davidson/diagonalize_CI_mono.irp.f +++ b/src/Davidson/diagonalize_CI_mono.irp.f @@ -16,7 +16,7 @@ if (diag_algorithm == "Davidson") then call davidson_diag(psi_det,CI_eigenvectors_mono,CI_electronic_energy, & - size(CI_eigenvectors_mono,1),N_det,N_states_diag,N_int,output_determinants) + size(CI_eigenvectors_mono,1),N_det,N_states,N_states_diag,N_int,output_determinants) else if (diag_algorithm == "Lapack") then diff --git a/src/Determinants/SC2.irp.f b/src/Determinants/SC2.irp.f deleted file mode 100644 index 4f321b87..00000000 --- a/src/Determinants/SC2.irp.f +++ /dev/null @@ -1,216 +0,0 @@ -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(output_determinants) - write(output_determinants,'(A)') '' - write(output_determinants,'(A)') 'CISD SC2' - write(output_determinants,'(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,Nint,output_determinants) - 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(output_determinants,'(A,I3)') 'SC2 Iteration ', iter - write(output_determinants,'(A)') '------------------' - write(output_determinants,'(A)') '' - write(output_determinants,'(A)') '===== ================' - write(output_determinants,'(A)') 'State Energy ' - write(output_determinants,'(A)') '===== ================' - do i=1,N_st - write(output_determinants,'(I5,1X,F16.10)') i, energies(i)+nuclear_repulsion - enddo - write(output_determinants,'(A)') '===== ================' - write(output_determinants,'(A)') '' - call write_double(output_determinants,(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(output_determinants) - deallocate (doubles,e_corr_array,H_jj_ref,H_jj_dressed, & - index_double, degree_exc, hij_double) - -end - - diff --git a/src/Determinants/options.irp.f b/src/Determinants/options.irp.f index d4283128..365021db 100644 --- a/src/Determinants/options.irp.f +++ b/src/Determinants/options.irp.f @@ -12,6 +12,9 @@ BEGIN_PROVIDER [ integer, N_states_diag ] else N_states_diag = N_states endif + if (N_states_diag < N_states) then + N_states_diag = N_states + endif call write_time(output_determinants) call write_int(output_determinants, N_states_diag, &