diff --git a/src/CISD/README.rst b/src/CISD/README.rst index 1a9c86d9..fbf1807f 100644 --- a/src/CISD/README.rst +++ b/src/CISD/README.rst @@ -23,12 +23,15 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`fill_h_apply_buffer_cisd `_ +`fill_h_apply_buffer_cisd `_ Fill the H_apply buffer with determiants for CISD -`h_apply_cisd `_ +`h_apply_cisd `_ Calls H_apply on the HF determinant and selects all connected single and double excitations (of the same symmetry). +`cisd `_ + Undocumented + diff --git a/src/CISD/SC2.irp.f b/src/CISD/SC2.irp.f new file mode 100644 index 00000000..330e3665 --- /dev/null +++ b/src/CISD/SC2.irp.f @@ -0,0 +1,421 @@ +subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint) + 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) + PROVIDE ref_bitmask_energy + 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,index_double(sze) + double precision :: hij_elec, e_corr_double,e_corr,diag_h_mat_elem,inv_c0 + double precision :: e_corr_array(sze),H_jj_ref(sze),H_jj_dressed(sze),hij_double(sze) + double precision :: e_corr_double_before,accu,cpu_2,cpu_1 + integer :: i_ok + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP SHARED(sze,N_st, & + !$OMP H_jj_ref,Nint,dets_in,u_in) & + !$OMP PRIVATE(i) + + !$OMP DO + 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) + if(degree==0)then + index_hf=i + else if (degree == 2)then + N_double += 1 + index_double(N_double) = 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) + index_double(N_double) = i + else if (degree == 1)then + call i_H_j(ref_bitmask,dets_in(1,1,i),Nint,hij_elec) + print*,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 + print*, 'E_corr = ',e_corr + print*, 'E_corr_double = ', e_corr_double + + converged = .False. + e_corr_double_before = e_corr_double + iter = 0 + do while (.not.converged) + + iter +=1 + print*,'SC2 iteration : ',iter + call cpu_time(cpu_1) + do i=1,sze + H_jj_dressed(i) = H_jj_ref(i) + if (i==index_hf)cycle + accu = 0.d0 + do j=1,N_double + call repeat_excitation(dets_in(1,1,i),ref_bitmask,dets_in(1,1,index_double(j)),i_ok,Nint) + if (i_ok==1)cycle! you check if the excitation is possible + accu += e_corr_array(j) + enddo + H_jj_dressed(i) += accu + enddo + + call cpu_time(cpu_2) + print*,'time for the excitations = ',cpu_2 - cpu_1 + print*,H_jj_ref(1),H_jj_ref(2) + print*,H_jj_dressed(1),H_jj_dressed(2) + print*,u_in(index_hf,1),u_in(index_double(1),1) + call davidson_diag_hjj(dets_in,u_in,H_jj_dressed,energies,dim_in,sze,N_st,Nint) + print*,u_in(index_hf,1),u_in(index_double(1),1) + 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 + print*,'E_corr = ',e_corr_double + print*,'delta E_corr =',e_corr_double - e_corr_double_before + converged = dabs(e_corr_double - e_corr_double_before) < 1.d-10 + if (converged) then + exit + endif + e_corr_double_before = e_corr_double + + enddo + + + +end + +subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint) + 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 + ! + ! 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(in) :: H_jj(dim_in) + 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) + + PROVIDE ref_bitmask_energy + + 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(W(1,k,iter),U(1,k,iter),H_jj,sze,dets_in,Nint) + 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 + + ! Diagonalize h + ! ------------- + call lapack_diag(lambda,y,h,N_st*davidson_sze_max,N_st*iter) + + ! Express eigenvectors of h in the determinant basis + ! -------------------------------------------------- + + ! call dgemm ( 'N','N', sze, N_st*iter, N_st, & + ! 1.d0, U(1,1,1), size(U,1), y(1,1,1,1), size(y,1)*size(y,2), & + ! 0.d0, U(1,1,iter+1), size(U,1) ) + 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) + enddo + + print '(I3,15(F16.8,x))', iter, lambda(1:N_st) + nuclear_repulsion + print '(3x,15(E16.5,x))', residual_norm(1:N_st) + + converged = maxval(residual_norm) < 1.d-10 + if (converged) then + exit + endif + + ! Davidson step + ! ------------- + + do k=1,N_st + do i=1,sze + U(i,k,iter+1) = 1.d0/(lambda(k) - H_jj(i)) * 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 + 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 + + deallocate ( & + kl_pairs, & + W, & + U, & + R, & + h, & + y, & + lambda & + ) +end + +subroutine repeat_excitation(key_in,key_1,key_2,i_ok,Nint) + use bitmasks + implicit none + integer(bit_kind), intent(in) :: key_in(Nint,2),key_1(Nint,2),key_2(Nint,2),Nint + integer,intent(out):: i_ok + integer :: ispin,i_hole,k_hole,j_hole,i_particl,k_particl,j_particl,i_trou,degree,exc(0:2,2,2) + double precision :: phase + i_ok = 1 + call get_excitation(key_1,key_2,exc,degree,phase,Nint) + integer :: h1,p1,h2,p2,s1,s2 + if(degree==2)then + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + ! first hole + k_hole = ishft(h1-1,-5)+1 + j_hole = h1-ishft(k_hole-1,5)-1 + if(iand(key_in(k_hole,s1),ibset(0,j_hole)).eq.0)then + i_ok = 0 + return + endif + + ! second hole + k_hole = ishft(h2-1,-5)+1 + j_hole = h2-ishft(k_hole-1,5)-1 + if(iand(key_in(k_hole,s2),ibset(0,j_hole)).eq.0)then + i_ok = 0 + return + endif + + ! first particle + k_particl = ishft(p1-1,-5)+1 + j_particl = p1-ishft(k_particl-1,5)-1 + if(iand(key_in(k_particl,s1),ibset(0,j_particl)).ne.0)then + i_ok = 0 + return + endif + + ! second particle + k_particl = ishft(p2-1,-5)+1 + j_particl = p2-ishft(k_particl-1,5)-1 + if(iand(key_in(k_particl,s2),ibset(0,j_particl)).ne.0)then + i_ok = 0 + return + endif + return + endif +end + diff --git a/src/CISD/cisd.irp.f b/src/CISD/cisd.irp.f index ba84b8c6..f3bb5488 100644 --- a/src/CISD/cisd.irp.f +++ b/src/CISD/cisd.irp.f @@ -6,6 +6,8 @@ program cisd double precision :: pt2(10), norm_pert(10), H_pert_diag + N_states = 3 + touch N_states call H_apply_cisd allocate(eigvalues(n_states),eigvectors(n_det,n_states)) print *, 'N_det = ', N_det @@ -20,7 +22,9 @@ program cisd print *, 'HF:', HF_energy print *, '---' do i = 1,1 - print *, 'energy(i) = ',eigvalues(i) + nuclear_repulsion + print *, 'energy(i) = ',eigvalues(i) + nuclear_repulsion + print *, 'E_corr = ',eigvalues(i) - ref_bitmask_energy enddo + call CISD_SC2(psi_det,psi_coef,eigvalues,size(psi_coef,1),N_det,N_states,N_int) deallocate(eigvalues,eigvectors) end diff --git a/src/CIS_dressed/CIS_providers.irp.f b/src/CIS_dressed/CIS_providers.irp.f new file mode 100644 index 00000000..b05f87cb --- /dev/null +++ b/src/CIS_dressed/CIS_providers.irp.f @@ -0,0 +1,334 @@ + + use bitmasks + BEGIN_PROVIDER [integer(bit_kind), psi_CIS,(N_int,2,size_psi_CIS)] + &BEGIN_PROVIDER [integer, psi_CIS_holes,(size_psi_CIS)] + &BEGIN_PROVIDER [integer, psi_CIS_particl,(size_psi_CIS)] + &BEGIN_PROVIDER [integer, psi_CIS_spin,(size_psi_CIS)] + &BEGIN_PROVIDER [integer, psi_CIS_adress,(n_core_cis+1:elec_alpha_num,elec_alpha_num+1:n_act_cis)] + &BEGIN_PROVIDER [double precision, H_CIS,(size_psi_CIS,size_psi_CIS)] + BEGIN_DOC + !key of the CIS-matrix + END_DOC + + + implicit none + integer :: a !control variable + integer :: i,j,k,l !variables for going over the occupied (i,j) and virutal (k,l) + integer :: key !key for CIS-matrix + integer :: i_hole,j_hole,ispin,l_particle,k_particle + double precision :: hij + + + do a=1,N_int + psi_CIS(a,1,1)=ref_bitmask(a,1) + psi_CIS(a,2,1)=ref_bitmask(a,2) + enddo + + psi_CIS_holes(1) = 0 + psi_CIS_particl(1) = 0 + psi_CIS_spin(1) = 0 + + !loop on particles: create a particle in k + do k=elec_alpha_num+1,n_act_cis + + !loop on holes: destroy a particle in i + do i=n_core_cis+1,elec_alpha_num + + !alpha spin + ispin=1 + + key=2*((k-elec_alpha_num-1)*(elec_alpha_num-n_core_cis) + i-n_core_cis) !index of such an excited determinant in the CIS WF + psi_CIS_adress(i,k)=key + + do a=1,N_int + psi_CIS(a,1,key)=ref_bitmask(a,1) + psi_CIS(a,2,key)=ref_bitmask(a,2) + enddo + + j_hole=ishft(i-1,-5)+1 + i_hole=i-ishft(j_hole-1,5)-1 + + psi_CIS(j_hole,ispin,key)=ibclr(psi_CIS(j_hole,ispin,key),i_hole) + + l_particle=ishft(k-1,-5)+1 + k_particle=k-ishft(l_particle-1,5)-1 + + psi_CIS(l_particle,ispin,key)=ibset(psi_CIS(l_particle,ispin,key),k_particle) + + psi_CIS_holes(key) = i + psi_CIS_particl(key) = k + psi_CIS_spin(key) = 1 + + !beta spin + ispin=2 + + key=key+1 + + do a=1,N_int + psi_CIS(a,1,key)=ref_bitmask(a,1) + psi_CIS(a,2,key)=ref_bitmask(a,2) + enddo + + j_hole=ishft(i-1,-5)+1 + i_hole=i-ishft(j_hole-1,5)-1 + + psi_CIS(j_hole,ispin,key)=ibclr(psi_CIS(j_hole,ispin,key),i_hole) + + l_particle=ishft(k-1,-5)+1 + k_particle=k-ishft(l_particle-1,5)-1 + psi_CIS_holes(key) = i + psi_CIS_particl(key) = k + psi_CIS_spin(key) = 2 + + psi_CIS(l_particle,ispin,key)=ibset(psi_CIS(l_particle,ispin,key),k_particle) + enddo + enddo + + !Building the CIS-matrix + double precision :: diag_H_mat_elem + + do key=1,size_psi_CIS + H_CIS(key,key)=diag_H_mat_elem(psi_CIS(1,1,key),N_int) + + do a=key+1,size_psi_CIS + call i_H_j(psi_CIS(1,1,a),psi_CIS(1,1,key),N_int,hij) + + H_CIS(key,a)=hij + H_CIS(a,key)=hij + enddo + enddo + + END_PROVIDER + + + BEGIN_PROVIDER[double precision, eigenvalues_CIS,(n_state_CIS)] + &BEGIN_PROVIDER[double precision, coefs_CIS, (size_psi_CIS,n_state_CIS)] + use bitmasks + + BEGIN_DOC + !the first states of the CIS matrix + END_DOC + + implicit none + integer :: i,j,k + double precision :: eigvalues(size_psi_CIS),eigvectors(size_psi_CIS,size_psi_CIS) + double precision :: coefs_tmp(size_psi_CIS) + double precision :: s2 + + + !Diagonalisation of CIS-matrix + call lapack_diag(eigvalues,eigvectors,H_CIS,size_psi_CIS,size_psi_CIS) + + do i = 1,n_state_CIS + eigenvalues_CIS(i) = eigvalues(i) + + do k=1,size_psi_CIS + + if (dabs(eigvectors(k,i)).ge.10.d-2) then + write(11,*),'k,i,eigenvectors(k,i)=',k,i,eigvectors(k,i) + write(11,*),'hole,particl,spin:',psi_CIS_holes(k),psi_CIS_particl(k),psi_CIS_spin(k) + write(11,*),'' + endif + coefs_tmp(k) = eigvectors(k,i) + coefs_CIS(k,i)=eigvectors(k,i) + enddo + call get_s2_u0(psi_CIS,coefs_tmp,size_psi_CIS,size_psi_CIS,s2) + enddo + + END_PROVIDER + + + BEGIN_PROVIDER [double precision, eigenvalues_CIS_dress_D,(n_state_CIS)] +&BEGIN_PROVIDER [double precision, s_2_CIS_dress_D,(n_state_CIS)] +&BEGIN_PROVIDER [double precision, eigenvectors_CIS_dress_D,(size_psi_CIS,n_state_CIS)] +&BEGIN_PROVIDER [double precision, overlap_D ] + use bitmasks + + BEGIN_DOC + !The first states of the CIS matrix dressed by the doubles + END_DOC + implicit none + double precision,allocatable :: delta_H_matrix_doub(:,:) + double precision,allocatable :: eigvalues(:),eigvectors(:,:) + double precision :: overlap,max_overlap,s2 + integer :: i_overlap,i,j,k + allocate (delta_H_matrix_doub(size_psi_CIS,size_psi_CIS)) + allocate(eigvalues(size_psi_CIS),eigvectors(size_psi_CIS,size_psi_CIS)) + do i = 1,n_state_CIS + call dress_by_doubles(eigenvalues_CIS(i),coefs_CIS(1,i),delta_H_matrix_doub,size_psi_CIS) !dressing of the Doubles + do j = 1,size_psi_CIS + do k = 1,size_psi_CIS + delta_H_matrix_doub(j,k) += H_CIS(j,k) + enddo + enddo + call lapack_diag(eigvalues,eigvectors,delta_H_matrix_doub,size_psi_CIS,size_psi_CIS) + + ! state following + max_overlap = 0.d0 + do k = 1, size_psi_CIS + overlap = 0.d0 + do j = 1,size_psi_CIS + overlap += eigvectors(j,k)*coefs_CIS(j,i) + enddo + if(dabs(overlap).gt.max_overlap)then + max_overlap = dabs(overlap) + i_overlap = k + endif + ! + enddo + ! print*,'overlap = ',max_overlap + overlap_D=max_overlap + do k = 1,size_psi_CIS + eigenvectors_CIS_dress_D(k,i) = eigvectors(k,i_overlap) + if (dabs(eigvectors(k,i_overlap)).ge.10.d-2) then + write(12,*),'k,i,eigenvectors(k,i)=',k,i,eigvectors(k,i_overlap) + write(12,*),'hole,particl,spin:',psi_CIS_holes(k),psi_CIS_particl(k),psi_CIS_spin(k) + write(12,*),'' + endif + enddo + call get_s2_u0(psi_CIS,eigenvectors_CIS_dress_D(1,i),size_psi_CIS,size_psi_CIS,s2) + s_2_CIS_dress_D(i) = s2 + eigenvalues_CIS_dress_D(i) = eigvalues(i_overlap) + enddo + + END_PROVIDER + + + + BEGIN_PROVIDER [double precision, eigenvalues_CIS_dress_D_dt,(n_state_CIS)] +&BEGIN_PROVIDER [double precision, s_2_CIS_dress_D_dt,(n_state_CIS)] +&BEGIN_PROVIDER [double precision, eigenvectors_CIS_dress_D_dt,(size_psi_CIS,n_state_CIS) +&BEGIN_PROVIDER [double precision, overlap_Ddt] + use bitmasks + + BEGIN_DOC + !The first states of the CIS matrix dressed by the doubles and the disconnected triples + END_DOC + implicit none + double precision,allocatable :: delta_H_matrix_doub(:,:) + double precision,allocatable :: eigvalues(:),eigvectors(:,:) + double precision :: overlap,max_overlap,s2 + integer :: i_overlap,i,j,k + allocate (delta_H_matrix_doub(size_psi_CIS,size_psi_CIS)) + allocate(eigvalues(size_psi_CIS),eigvectors(size_psi_CIS,size_psi_CIS)) + do i = 1,n_state_CIS + call dress_by_doubles(eigenvalues_CIS(i),coefs_CIS(1,i),delta_H_matrix_doub,size_psi_CIS) !dressing of the Doubles + + do j = 1,size_psi_CIS + do k = 1,size_psi_CIS + delta_H_matrix_doub(j,k) += H_CIS(j,k) + enddo + delta_H_matrix_doub(j,j) += dress_T_discon_array_CIS(j) + enddo + call lapack_diag(eigvalues,eigvectors,delta_H_matrix_doub,size_psi_CIS,size_psi_CIS) + + ! state following + max_overlap = 0.d0 + do k = 1, size_psi_CIS + overlap = 0.d0 + do j = 1,size_psi_CIS + overlap += eigvectors(j,k)*coefs_CIS(j,i) + enddo + if(dabs(overlap).gt.max_overlap)then + max_overlap = dabs(overlap) + i_overlap = k + endif + ! + enddo + ! print*,'overlap = ',max_overlap + overlap_Ddt=max_overlap + do k = 1,size_psi_CIS + eigenvectors_CIS_dress_D_dt(k,i) = eigvectors(k,i_overlap) + if (dabs(eigvectors(k,i_overlap)).ge.10.d-2) then + write(13,*),'k,i,eigenvectors(k,i)=',k,i,eigvectors(k,i_overlap) + write(13,*),'hole,particl,spin:',psi_CIS_holes(k),psi_CIS_particl(k),psi_CIS_spin(k) + write(13,*),'' + endif + enddo + call get_s2_u0(psi_CIS,eigenvectors_CIS_dress_D_dt(1,i),size_psi_CIS,size_psi_CIS,s2) + s_2_CIS_dress_D_dt(i) = s2 + eigenvalues_CIS_dress_D_dt(i) = eigvalues(i_overlap) + enddo + + END_PROVIDER + +!BEGIN_PROVIDER [double precision, eigenvalues_CIS_dress_tot,(n_state_CIS)] +!BEGIN_PROVIDER [double precision, s_2_CIS_dress_tot,(n_state_CIS)] +!BEGIN_PROVIDER [double precision, eigenvectors_CIS_dress_tot,(size_psi_CIS,n_state_CIS)] +!BEGIN_PROVIDER [double precision, overlap_tot] + +!BEGIN_DOC +!!The first states of the CIS matrix dressed by the doubles +!END_DOC +!implicit none +!double precision,allocatable :: delta_H_matrix_doub(:,:) +!double precision,allocatable :: eigvalues(:),eigvectors(:,:) +!double precision,allocatable :: delta_H_trip(:,:) +!double precision :: overlap,max_overlap,s2,average_eigvalue +!integer :: i_overlap,i,j,k,m +!allocate (delta_H_matrix_doub(size_psi_CIS,size_psi_CIS),delta_H_trip(size_psi_CIS,size_psi_CIS) ) +!allocate(eigvalues(size_psi_CIS),eigvectors(size_psi_CIS,size_psi_CIS)) +! do i = 1,n_state_CIS +! call dress_by_doubles(eigenvalues_CIS(i),coefs_CIS(1,i),delta_H_matrix_doub,size_psi_CIS) !dressing of the Doubles +! call dress_T_con(eigenvalues_CIS(i),delta_H_trip,size_psi_CIS) +! do j = 1,size_psi_CIS +! do k = 1,size_psi_CIS +! delta_H_matrix_doub(j,k) += H_CIS(j,k) +! delta_H_matrix_doub(j,k) += delta_H_trip(j,k) +! enddo +! delta_H_matrix_doub(j,j) += dress_T_discon_array_CIS(j) +! enddo +! call lapack_diag(eigvalues,eigvectors,delta_H_matrix_doub,size_psi_CIS,size_psi_CIS) +! do m=1,n_state_CIS +! write(12,*),'m,eigvalues(m)',m,eigvalues(m) +! enddo +! ! state following +! max_overlap = 0.d0 +! do k = 1, size_psi_CIS +! overlap = 0.d0 +! do j = 1,size_psi_CIS +! overlap += eigvectors(j,k)*coefs_CIS(j,i) +! enddo +! if(dabs(overlap).gt.max_overlap)then +! max_overlap = dabs(overlap) +! i_overlap = k +! endif +! ! +! enddo +!! print*,'overlap = ',max_overlap +! overlap_tot=max_overlap +! do k = 1,size_psi_CIS +! eigenvectors_CIS_dress_tot(k,i) = eigvectors(k,i_overlap) +! enddo +! call get_s2_u0(psi_CIS,eigenvectors_CIS_dress_tot(1,i),size_psi_CIS,size_psi_CIS,s2) +! s_2_CIS_dress_tot(i) = s2 +! eigenvalues_CIS_dress_tot(i) = eigvalues(i_overlap) +! enddo + +!END_PROVIDER + + + + + + + + + BEGIN_PROVIDER [double precision, diag_elements, (size_psi_CIS)] + use bitmasks + + + BEGIN_DOC + !Array of the energy of the CIS determinants ordered in the CIS matrix + END_DOC + + implicit none + double precision :: hij + integer :: i + + do i = 1, size_psi_CIS + call i_H_j(psi_CIS(1,1,i),psi_CIS(1,1,i),N_int,hij) + diag_elements(i) = hij + enddo + + END_PROVIDER diff --git a/src/CIS_dressed/MP2.irp.f b/src/CIS_dressed/MP2.irp.f new file mode 100644 index 00000000..cfcb3369 --- /dev/null +++ b/src/CIS_dressed/MP2.irp.f @@ -0,0 +1,1611 @@ + BEGIN_PROVIDER [double precision, MP2_corr_energy] + &BEGIN_PROVIDER [double precision, p_imp,(elec_alpha_num+1:n_act_cis)] + &BEGIN_PROVIDER [double precision, h_imp,(n_core_cis+1:elec_alpha_num)] + &BEGIN_PROVIDER [double precision, hp_imp,(n_core_cis+1:elec_alpha_num,elec_alpha_num+1:n_act_cis)] + + BEGIN_DOC + !Calculation of the MP2 correlation energy (MP2_corr_energy) + !and calculation of the contribution of the disconnected Triples on the + !Singles, via the impossible (p_imp, h_imp, hp_imp) + END_DOC + + implicit none + integer :: i,j,k,l !variables for going over the occupied (i,j) and virutal (k,l) + double precision :: direct,exchg,hij !calculating direct, exchange and total contribution + double precision :: get_mo_bielec_integral + double precision :: e_i,e_k,e_j,e_l !epsilons of i,j,k,l + double precision :: delta_e_ik,delta_e_ikj + double precision :: delta_e !delta epsilons + + + + MP2_corr_energy=0.d0 + + do k =elec_beta_num + 1, n_act_cis + p_imp(k)=0.d0 + enddo + + + do i=n_core_cis+1,elec_alpha_num + h_imp(i)=0.d0 + + e_i=diagonal_Fock_matrix_mo(i) + + do k=elec_alpha_num+1,n_act_cis + hp_imp(i,k)=0.d0 + + e_k=diagonal_Fock_matrix_mo(k) + delta_e_ik=e_i-e_k + + !same spin contribution for MP2_corr_energy + do j=i+1,elec_alpha_num + e_j=diagonal_Fock_matrix_mo(j) + delta_e_ikj=delta_e_ik+e_j + + !same spin contribution for MP2_corr_energy and a part of the contribution to p_imp and h_imp + do l=k+1,n_act_cis + e_l=diagonal_Fock_matrix_mo(l) + delta_e=delta_e_ikj-e_l + + direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + exchg =get_mo_bielec_integral(i,j,l,k,mo_integrals_map) + + hij=(direct-exchg)*(direct-exchg)/delta_e + + MP2_corr_energy=MP2_corr_energy+2*hij + p_imp(k)=p_imp(k)+hij + h_imp(i)=h_imp(i)+hij + + enddo + + !remaining same spin contribution for p_imp + do l=elec_alpha_num+1,k-1 + e_l=diagonal_Fock_matrix_mo(l) + delta_e=delta_e_ikj-e_l + + direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + exchg =get_mo_bielec_integral(i,j,l,k,mo_integrals_map) + + hij=(direct-exchg)*(direct-exchg)/delta_e + + p_imp(k)=p_imp(k)+hij + + enddo + enddo + + !remaining same spin contribution for h_imp + do j=n_core_cis+1,i-1 + e_j=diagonal_Fock_matrix_mo(j) + delta_e_ikj=delta_e_ik+e_j + + do l=k+1,n_act_cis + e_l=diagonal_Fock_matrix_mo(l) + delta_e=delta_e_ikj-e_l + + direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + exchg =get_mo_bielec_integral(i,j,l,k,mo_integrals_map) + + hij=(direct-exchg)*(direct-exchg)/delta_e + + h_imp(i)=h_imp(i)+hij + + enddo + enddo + + !same spin contribution for hp_imp + do j=n_core_cis+1,elec_alpha_num + e_j=diagonal_Fock_matrix_mo(j) + delta_e_ikj=delta_e_ik+e_j + + do l=elec_alpha_num+1,n_act_cis + e_l=diagonal_Fock_matrix_mo(l) + delta_e=delta_e_ikj-e_l + + direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + exchg =get_mo_bielec_integral(i,j,l,k,mo_integrals_map) + + hij=(direct-exchg)*(direct-exchg)/delta_e + + hp_imp(i,k)=hp_imp(i,k)+hij + + enddo + enddo + + !different spin contribution + do j=n_core_cis+1,elec_beta_num + e_j=diagonal_Fock_matrix_mo(j) + delta_e_ikj=delta_e_ik+e_j + + do l=elec_beta_num+1,n_act_cis + e_l=diagonal_Fock_matrix_mo(l) + delta_e=delta_e_ikj-e_l + + direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + + hij=direct*direct/delta_e + + MP2_corr_energy=MP2_corr_energy+hij + p_imp(k)=p_imp(k)+hij + h_imp(i)=h_imp(i)+hij + hp_imp(i,k)=hp_imp(i,k)+hij + + enddo + enddo + enddo + enddo + + print*,'MP2 correlation energy=',MP2_corr_energy + + END_PROVIDER + + + BEGIN_PROVIDER [double precision, EN2_corr_energy] + &BEGIN_PROVIDER [double precision, p_imp_EN,(elec_alpha_num+1:n_act_cis)] + &BEGIN_PROVIDER [double precision, h_imp_EN,(n_core_cis+1:elec_alpha_num)] + &BEGIN_PROVIDER [double precision, hp_imp_EN,(n_core_cis+1:elec_alpha_num,elec_alpha_num+1:n_act_cis)] + + BEGIN_DOC + !Calculation of the EN2 correlation energy (EN2_corr_energy) + !and calculation of the contribution of the disconnected Triples on the + !Singles, via the impossible (p_imp_EN, h_imp_EN, hp_imp_EN) + END_DOC + + implicit none + integer :: i,j,k,l !variables for going over the occupied (i,j) and virutal (k,l) + double precision :: direct,exchg,hij !calculating direct, exchange and total contribution + double precision :: get_mo_bielec_integral + double precision :: e_i,e_k,e_j,e_l !epsilons of i,j,k,l + double precision :: delta_e_ik,delta_e_ikj + double precision :: delta_e !delta epsilons + double precision :: delta_e_tmp,H_jj_total + integer, allocatable :: key_in(:,:) + integer, allocatable :: key_out(:,:) + integer :: ispin1,ispin2,i_ok + allocate(key_in(N_int,2)) + allocate(key_out(N_int,2)) + + print*,'EN2_corr_energy' + + EN2_corr_energy=0.d0 + +do i=n_core_cis+1,elec_alpha_num + h_imp_EN(i)=0.d0 + do k =elec_beta_num + 1, n_act_cis + p_imp_EN(k)=0.d0 + hp_imp_EN(i,k)=0.d0 + enddo +enddo + + do i=n_core_cis+1,elec_alpha_num +! h_imp_EN(i)=0.d0 + + e_i=diagonal_Fock_matrix_mo(i) + + do k=elec_alpha_num+1,n_act_cis +! hp_imp_EN(i,k)=0.d0 + + e_k=diagonal_Fock_matrix_mo(k) + delta_e_ik=e_i-e_k + + !same spin contribution for EN2_corr_energy + do j=i+1,elec_alpha_num + e_j=diagonal_Fock_matrix_mo(j) + delta_e_ikj=delta_e_ik+e_j + + !same spin contribution for EN2_corr_energy and a part of the contribution to p_imp_EN,h_imp_EN + do l=k+1,n_act_cis + e_l=diagonal_Fock_matrix_mo(l) + delta_e=delta_e_ikj-e_l + delta_e=delta_e-mo_bielec_integral_jj_anti(i,j) - & + mo_bielec_integral_jj_anti(k,l) + delta_e=delta_e+mo_bielec_integral_jj_anti(i,k) + & + mo_bielec_integral_jj_anti(i,l) + & + mo_bielec_integral_jj_anti(j,k) + & + mo_bielec_integral_jj_anti(j,l) + ! ispin1 = 1 + ! ispin2 = 1 + ! call diexcitation(i,j,k,l,ispin1,ispin2,ref_bitmask,key_out,i_ok,N_int) + ! delta_e_tmp = HF_energy - H_jj_total(key_out) + ! if(dabs(delta_e_tmp-delta_e).gt.1.d-10)then + ! print*,'same spin first' + ! print*,'delta_e_SSS = ',delta_e_tmp-delta_e + ! stop + ! endif + ! + direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + exchg=get_mo_bielec_integral(i,j,l,k,mo_integrals_map) + + hij=(direct-exchg)*(direct-exchg) + hij=0.5d0*(-delta_e-dsqrt(delta_e*delta_e+4.d0*hij)) + ! ispin1 = 1 + ! ispin2 = 1 + ! call diexcitation(i,j,k,l,ispin1,ispin2,ref_bitmask,key_out,i_ok,N_int) + ! delta_e_tmp = HF_energy - H_jj_total(key_out) + + EN2_corr_energy=EN2_corr_energy+2*hij + p_imp_EN(k)=p_imp_EN(k)+hij + h_imp_EN(i)=h_imp_EN(i)+hij + + p_imp_EN(l)=p_imp_EN(l)+hij + h_imp_EN(j)=h_imp_EN(j)+hij + + enddo +enddo +! !remaining same spin contribution for p_imp_EN + +! do l=elec_alpha_num+1,k-1 +! e_l=diagonal_Fock_matrix_mo(l) +! delta_e=delta_e_ikj-e_l +! delta_e=e_i + e_j - e_k -e_l +! delta_e=delta_e-mo_bielec_integral_jj_anti(i,j) - & +! mo_bielec_integral_jj_anti(k,l) +! delta_e=delta_e+mo_bielec_integral_jj_anti(i,k) + & +! mo_bielec_integral_jj_anti(i,l) + & +! mo_bielec_integral_jj_anti(j,k) + & +! mo_bielec_integral_jj_anti(j,l) +! e_l=diagonal_Fock_matrix_mo(l) +! delta_e=delta_e_ikj-e_l +! delta_e=delta_e-mo_bielec_integral_jj_anti(i,j) - & +! mo_bielec_integral_jj_anti(k,l) +! delta_e=delta_e+mo_bielec_integral_jj_anti(i,k) + & +! mo_bielec_integral_jj_anti(i,l) + & +! mo_bielec_integral_jj_anti(j,k) + & +! mo_bielec_integral_jj_anti(j,l) + + + + + +!! ispin1 = 1 +!! ispin2 = 1 +!! call diexcitation(i,j,k,l,ispin1,ispin2,ref_bitmask,key_out,i_ok,N_int) +!! delta_e_tmp = HF_energy - H_jj_total(key_out) +!! if(dabs(delta_e_tmp-delta_e).gt.1.d-10)then +!! call print_key(key_out) +!! call print_key(ref_bitmask) +!! print*,i,k,j,l +!! print*,'same spin middle' +!! print*,'delta_e_SSS = ',delta_e_tmp-delta_e +!! print*,'delta_e_SSS = ',delta_e_tmp,delta_e +!! stop +!! endif + +! direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) +! exchg=get_mo_bielec_integral(i,j,l,k,mo_integrals_map) + +! hij=(direct-exchg)*(direct-exchg) +! hij=0.5d0*(-delta_e-dsqrt(delta_e*delta_e+4.d0*hij)) + +! p_imp_EN(k)=p_imp_EN(k)+hij + +! enddo +! enddo + +! !remaining same spin contribution for h_imp_EN +! do j=n_core_cis+1,i-1 +! e_j=diagonal_Fock_matrix_mo(j) +! delta_e_ikj=delta_e_ik+e_j + +! do l=k+1,n_act_cis +! e_l=diagonal_Fock_matrix_mo(l) +! delta_e=delta_e_ikj-e_l +! delta_e=delta_e-mo_bielec_integral_jj_anti(i,j) - & +! mo_bielec_integral_jj_anti(k,l) +! delta_e=delta_e+mo_bielec_integral_jj_anti(i,k) + & +! mo_bielec_integral_jj_anti(i,l) + & +! mo_bielec_integral_jj_anti(j,k) + & +! mo_bielec_integral_jj_anti(j,l) +!! ispin2 = 1 +!! call diexcitation(i,j,k,l,ispin1,ispin2,ref_bitmask,key_out,i_ok,N_int) +!! delta_e_tmp = HF_energy - H_jj_total(key_out) +!! if(dabs(delta_e_tmp-delta_e).gt.1.d-10)then +!! print*,'same spin third ' +!! print*,'delta_e_SSS = ',delta_e_tmp-delta_e +!! stop +!! endif + +! direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) +! exchg=get_mo_bielec_integral(i,j,l,k,mo_integrals_map) + +! hij=(direct-exchg)*(direct-exchg) +! hij=0.5d0*(-delta_e-dsqrt(delta_e*delta_e+4.d0*hij)) + +! h_imp_EN(i)=h_imp_EN(i)+hij + +! enddo +! enddo + + !same spin contribution for hp_imp_EN + + do j=n_core_cis+1,elec_alpha_num + if(j==i)cycle + e_j=diagonal_Fock_matrix_mo(j) + delta_e_ikj=delta_e_ik+e_j + + do l=elec_alpha_num+1,n_act_cis + if(l==k)cycle + e_l=diagonal_Fock_matrix_mo(l) + delta_e=delta_e_ikj-e_l + delta_e=delta_e-mo_bielec_integral_jj_anti(i,j) - & + mo_bielec_integral_jj_anti(k,l) + delta_e=delta_e+mo_bielec_integral_jj_anti(i,k) + & + mo_bielec_integral_jj_anti(i,l) + & + mo_bielec_integral_jj_anti(j,k) + & + mo_bielec_integral_jj_anti(j,l) +!!! ispin1 = 1 +!!! ispin2 = 1 +!!! call diexcitation(i,j,k,l,ispin1,ispin2,ref_bitmask,key_out,i_ok,N_int) +!!! delta_e_tmp = HF_energy - H_jj_total(key_out) +!!! if(dabs(delta_e_tmp-delta_e).gt.1.d-10)then +!!! print*,'same spin fourth' +!!! print*,'delta_e_SSS = ',delta_e_tmp-delta_e +!!! call print_key(key_out) +!!! call print_key(ref_bitmask) +!!! print*,i,k,j,l +!!! stop +!!! endif + + direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + exchg=get_mo_bielec_integral(i,j,l,k,mo_integrals_map) + + hij=(direct-exchg)*(direct-exchg) + hij = 0.5d0 * (-delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij)) + + hp_imp_EN(i,k)=hp_imp_EN(i,k)+hij + + enddo + enddo + + !different spin contribution + do j=n_core_cis+1,elec_beta_num + e_j=diagonal_Fock_matrix_mo(j) + delta_e_ikj=delta_e_ik+e_j + + do l=elec_beta_num+1,n_act_cis + e_l=diagonal_Fock_matrix_mo(l) + delta_e=delta_e_ikj-e_l + delta_e=delta_e-mo_bielec_integral_jj(i,j) - & + mo_bielec_integral_jj(k,l) + delta_e=delta_e+mo_bielec_integral_jj_anti(i,k) + & + mo_bielec_integral_jj_anti(j,l) + & + mo_bielec_integral_jj(i,l) + & + mo_bielec_integral_jj(j,k) +!! ispin1 = 2 +!! ispin2 = 1 +!! call diexcitation(i,j,k,l,ispin1,ispin2,ref_bitmask,key_out,i_ok,N_int) +!! delta_e_tmp = HF_energy - H_jj_total(key_out) +!! if(dabs(delta_e_tmp-delta_e).gt.1.d-10)then +!! print*,'different spin ' +!! print*,'delta_e_SSS = ',delta_e_tmp-delta_e +!! stop +!! endif + + + direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + + hij=direct*direct + hij=0.5d0*(-delta_e-dsqrt(delta_e*delta_e+4.d0*hij)) + + EN2_corr_energy=EN2_corr_energy+hij + p_imp_EN(k)=p_imp_EN(k)+hij + h_imp_EN(i)=h_imp_EN(i)+hij + hp_imp_EN(i,k)=hp_imp_EN(i,k)+hij + + enddo + enddo + enddo + enddo + + print*,'EN correlation energy=',EN2_corr_energy + print*,'EN correlation energy=',EN2_corr_energy + + END_PROVIDER + + + + + BEGIN_PROVIDER [integer, size_psi_CIS] + + BEGIN_DOC + !Definition of the size of the CIS vector + END_DOC + + implicit none + size_psi_CIS=(elec_alpha_num-n_core_cis)*(n_act_cis-elec_alpha_num)*2+1 !!!Is this still correct for open shell systems??? + print*,'size_psi_CIS = ',size_psi_CIS + + END_PROVIDER + + + BEGIN_PROVIDER [double precision, diag_elements_sorted, (size_psi_CIS)] + &BEGIN_PROVIDER [integer , diag_elements_sorted_index, (size_psi_CIS)] + BEGIN_DOC + !Array of the energy of the CIS determinants sorted by energy and + !Index in the CIS matrix + END_DOC + + implicit none + integer :: iorder(size_psi_CIS),i + + print*,'diag_elements_sorted' + + do i = 1,size_psi_CIS + diag_elements_sorted(i) = diag_elements(i) + iorder(i) = i + enddo + + call dsort(diag_elements_sorted,iorder,size_psi_CIS) + + do i = 1,size_psi_CIS + diag_elements_sorted(i) = diag_elements(iorder(i)) + diag_elements_sorted_index(i) = iorder(i) + enddo + + END_PROVIDER + + + + + + + BEGIN_PROVIDER [double precision, eigenvalues_dressed_CIS,(n_state_CIS)] + + BEGIN_DOC + !The first states of the dressed CIS matrix + END_DOC + + implicit none + integer :: i,j,k + double precision :: s2,e_corr,tmp + double precision,allocatable :: coefs_tmp(:) + double precision,allocatable :: eigvalues(:),eigvectors(:,:) + double precision,allocatable :: delta_H_trip(:,:) + double precision,allocatable :: delta_H_matrix_doub(:,:) + double precision,allocatable :: delta_H_matrix_dpdiscon(:,:) + double precision,allocatable :: delta_H_matrix(:,:) + double precision :: eigenvalues_dressed_CIS_D(n_state_CIS) + double precision :: eigenvalues_dressed_CIS_DT(n_state_CIS) + double precision,allocatable :: eigvald(:),eigvecd(:,:) + allocate(eigvalues(size_psi_CIS),eigvectors(size_psi_CIS,size_psi_CIS),delta_H_trip(size_psi_CIS,size_psi_CIS),delta_H_matrix_doub(size_psi_CIS,size_psi_CIS),delta_H_matrix(size_psi_CIS,size_psi_CIS)) + allocate (eigvald(size_psi_CIS),eigvecd(size_psi_CIS,size_psi_CIS),coefs_tmp(size_psi_CIS),delta_H_matrix_dpdiscon(size_psi_CIS,size_psi_CIS) ) + + print*,'eigenvalues_dressed_CIS' + + provide coefs_CIS + do i = 1,n_state_CIS + + write(66,*),'i=',i + + call dress_by_doubles(eigenvalues_CIS(i),coefs_CIS(1,i),delta_H_matrix,size_psi_CIS) !dressing of the Doubles + do j = 1,size_psi_CIS + do k = 1,size_psi_CIS + delta_H_matrix(j,k)=delta_H_matrix(j,k)+H_CIS(j,k) + + delta_H_matrix_doub(j,k)=delta_H_matrix(j,k) + delta_H_matrix(j,k)=delta_H_matrix(j,k)+delta_H_trip(j,k) + delta_H_matrix_dpdiscon(j,k)=delta_H_matrix(j,k) + enddo + delta_H_matrix(j,j)=delta_H_matrix(j,j)+dress_T_discon_array_CIS(j) !dressing by the disconnected Triples + enddo + + + + + call lapack_diag(eigvalues,eigvectors,delta_H_matrix_doub,size_psi_CIS,size_psi_CIS) + double precision :: overlap + double precision :: max_overlap + integer :: i_overlap + max_overlap = 0.d0 + do k = 1, size_psi_CIS + overlap = 0.d0 + do j = 1,size_psi_CIS + overlap += eigvectors(j,k)*coefs_CIS(j,i) + enddo + if(dabs(overlap).gt.max_overlap)then + max_overlap = dabs(overlap) + i_overlap = k + endif + ! + enddo +! print*,'i_overlap = ',i_overlap + eigenvalues_dressed_CIS_D(i) = eigvalues(i_overlap) + + do j = 1,size_psi_CIS + coefs_tmp(j) = eigvectors(j,i_overlap) + + enddo + call get_s2_u0(psi_CIS,coefs_tmp,size_psi_CIS,size_psi_CIS,s2) +! print*,'s2(D)=',s2 +! print*,'' +! print*,'diagonal elements sorted(i) =',diag_elements_sorted(i) +! print*,'eigenvalues of the CIS(i) =',eigenvalues_CIS(i) +! print*,'eigenvalues dressed by Doubles (D) =',eigenvalues_dressed_CIS_D(i) + +! write(12,*),'' +! write(12,*),'i=',i +! write(12,*),'eigenvalues of the CIS(i) =',eigenvalues_CIS(i) +! write(12,*),'diagonal elements sorted(i)=',diag_elements_sorted(i) +! write(12,*),'' +! write(12,*),'s2(D)=',s2 +! write(12,*),'eigenvalues dressed by Doubles (D) =',eigenvalues_dressed_CIS_D(i) + + write(66,*),'CIS-HF ',eigenvalues_CIS(i)-eigenvalues_CIS(1) + write(66,*),'Doub-HF ',eigenvalues_dressed_CIS_D(i)-eigenvalues_dressed_CIS_D(1) + + call lapack_diag(eigvalues,eigvectors,delta_H_matrix_dpdiscon,size_psi_CIS,size_psi_CIS) + max_overlap = 0.d0 + do k = 1, size_psi_CIS + overlap = 0.d0 + do j = 1,size_psi_CIS + overlap += eigvectors(j,k)*coefs_CIS(j,i) + enddo + if(dabs(overlap).gt.max_overlap)then + max_overlap = dabs(overlap) + i_overlap = k + endif + ! + enddo +! print*,'i_overlap = ',i_overlap + + + eigenvalues_dressed_CIS_DT(i) = eigvalues(i_overlap) + + do j = 1,size_psi_CIS + coefs_tmp(j) = eigvectors(j,i_overlap) + enddo + + call get_s2_u0(psi_CIS,coefs_tmp,size_psi_CIS,size_psi_CIS,s2) + +! write(12,*),'s2(D+cT)=',s2 +! write(12,*),'eigenvalues dressed by D and connected Triples=',eigenvalues_dressed_CIS_DT(i) +! print*,'s2(DcT)=',s2 +! print*,'eigenvalues dressed D and cT=',eigenvalues_dressed_CIS_DT(i) + + write(66,*),'D+cT-HF ',eigenvalues_dressed_CIS_DT(i)-eigenvalues_dressed_CIS_DT(1) + + + + call lapack_diag(eigvalues,eigvectors,delta_H_matrix,size_psi_CIS,size_psi_CIS) + do k = 1, size_psi_CIS + overlap = 0.d0 + do j = 1,size_psi_CIS + overlap += eigvectors(j,k)*coefs_CIS(j,i) + enddo + if(dabs(overlap).gt.max_overlap)then + max_overlap = dabs(overlap) + i_overlap = k + endif + ! + enddo +! print*,'i_overlap = ',i_overlap + + + eigenvalues_dressed_CIS(i) = eigvalues(i_overlap) + + + do j = 1,size_psi_CIS + coefs_tmp(j) = eigvectors(j,i_overlap) + enddo + call get_s2_u0(psi_CIS,coefs_tmp,size_psi_CIS,size_psi_CIS,s2) + +! print*,'s2(Tot)=',s2 +! print*,'eigenvalues dressed =',eigenvalues_dressed_CIS(i) + +! write(12,*),'' +! write(12,*),'s2(tot)=',s2 +! write(12,*),'eigenvalues of dressed CIS(i)=',eigenvalues_dressed_CIS(i) + + write(66,*)'Eigval-HF ',eigenvalues_dressed_CIS(i)-eigenvalues_dressed_CIS(1) + enddo + + deallocate(eigvalues,eigvectors,delta_H_trip,delta_H_matrix_doub,delta_H_matrix_dpdiscon) + deallocate (eigvald,eigvecd,coefs_tmp) +!call system ( "mkdir results_$(date +%d%m%y_%I%M)" ) +!call system ( "cp ../newcode/MP2.irp.f results_$(date +%d%m%y_%I%M)" ) +!call system ( "mv fort.10 Ecorr.out" ) +!call system ( "mv fort.11 CIS.out" ) +!call system ( "mv fort.12 Dressing.out" ) +!call system ( "mv fort.66 delta_e_states" ) +!call system ( "cp Ecorr.out results_$(date +%d%m%y_%I%M)" ) +!call system ( "cp CIS.out results_$(date +%d%m%y_%I%M)" ) +!call system ( "cp Dressing.out results_$(date +%d%m%y_%I%M)" ) +!call system ( "cp delta_e_states results_$(date +%d%m%y_%I%M)" ) +!call system ( "rm delta_e_states" ) +!call system ( "rm Ecorr.out" ) +!call system ( "rm CIS.out" ) +!call system ( "rm Dressing.out" ) + + END_PROVIDER + + + BEGIN_PROVIDER [double precision, dress_T_discon,(n_core_cis+1:elec_alpha_num,elec_alpha_num+1:n_act_cis)] + &BEGIN_PROVIDER [double precision, dress_T_discon_array_CIS,(size_psi_CIS)] + + BEGIN_DOC + !Calculation of the dressing by the disconnected Triples, via the impossible + END_DOC + + implicit none + integer :: i,k !variables for going over the occupied (i) and virutal (k) + integer :: key !key for CIS-matrix + + print*,'dress_T_discon' + + if(mp2_dressing)then + dress_T_discon_array_CIS(1)=MP2_corr_energy + + do i=n_core_cis+1,elec_alpha_num + do k=elec_alpha_num+1,n_act_cis + key=psi_CIS_adress(i,k) + + dress_T_discon(i,k)=MP2_corr_energy-p_imp(k)-h_imp(i)+hp_imp(i,k) + + dress_T_discon_array_CIS(key) = dress_T_discon(i,k) + dress_T_discon_array_CIS(key+1) = dress_T_discon(i,k) + + enddo + enddo + + else !EN Dressing + dress_T_discon_array_CIS(1)=EN2_corr_energy + + do i=n_core_cis+1,elec_alpha_num + do k=elec_alpha_num+1,n_act_cis + key=psi_CIS_adress(i,k) + + dress_T_discon(i,k)=EN2_corr_energy-p_imp_EN(k)-h_imp_EN(i)+hp_imp_EN(i,k) + +!print*,'i,k,key',i,k,key +!print*,'EN2_corr_energy-p_imp_EN(k)-h_imp_EN(i)+hp_imp_EN(i,k)',EN2_corr_energy,p_imp_EN(k),h_imp_EN(i),hp_imp_EN(i,k) +!print*,'dress_T_discon',dress_T_discon(i,k) + + dress_T_discon_array_CIS(key) = dress_T_discon(i,k) + dress_T_discon_array_CIS(key+1) = dress_T_discon(i,k) + + enddo + enddo + end if + + END_PROVIDER + + + subroutine dress_by_doubles(ref_energy,coefs,delta_H_matrix,ndet) + + use bitmasks + implicit none + integer, intent(in) :: ndet + double precision, intent(in) :: ref_energy + double precision, intent(in) ::coefs(size_psi_CIS) + double precision, intent(out) :: delta_H_matrix(ndet,ndet) + + integer :: i,j,k,l !variables for going over the occupied (i,j) and virutal (k,l) + integer :: key !key for CIS-matrix + integer :: a,b,i_ok !control variables + integer(bit_kind) :: key_out(N_int,2) !Doubles + integer :: ispin1,ispin2 !spin variables + integer :: degree,i_count_connected + double precision :: hmD,hij,hij_array(ndet),hij_index(ndet) + double precision :: delta_hij + double precision :: delta_e , e_double!delta epsilons + double precision :: accu_H_mD,diag_H_mat_elem + double precision :: epsilon_D,s2 + + print*,'dress_by_doubles' + + double precision :: accu,accu_2 + accu = 0.d0 + accu_2 = 0.d0 + do i =1,ndet + accu += coefs(i) * coefs(i) + enddo + if (dabs(accu-1.d0) > 1.d-10)then + print*,'coefficients not normalized !!' + print*,accu + !stop + endif + + accu = 0.d0 + if(standard_doubles)then + delta_H_matrix = 0.d0 + + do i=n_core_cis+1,elec_alpha_num + + do k=elec_alpha_num+1,n_act_cis + + !same spin contribution + do j=i+1,elec_alpha_num + + do l=k+1,n_act_cis + !Alpha + ispin1=1 + ispin2=1 + + call diexcitation(i,j,k,l,ispin1,ispin2,ref_bitmask,key_out,i_ok,N_int) + + e_double =diag_H_mat_elem(key_out,N_int) + + delta_e = 1.d0 / (ref_energy - e_double) + + i_count_connected = 0 + + do a=2, ndet + call get_excitation_degree(key_out,psi_CIS(1,1,a),degree,N_int) + + if (degree /= 1 .and.degree /= 2) cycle + + call i_H_j(key_out,psi_CIS(1,1,a),N_int,hij) + + if (dabs(hij).le.1.d-10) cycle + + i_count_connected = i_count_connected +1 + + hij_array(i_count_connected) = hij + hij_index(i_count_connected) = a + + enddo + + do a=1,i_count_connected + delta_hij = hij_array(a)*hij_array(a) * delta_e + + delta_H_matrix(hij_index(a),hij_index(a))=delta_H_matrix(hij_index(a),hij_index(a))+delta_hij + + do b=a+1,i_count_connected + delta_hij = hij_array(a)*hij_array(b) * delta_e + + delta_H_matrix(hij_index(a),hij_index(b))=delta_H_matrix(hij_index(a),hij_index(b))+delta_hij + delta_H_matrix(hij_index(b),hij_index(a))=delta_H_matrix(hij_index(b),hij_index(a))+delta_hij + + enddo + enddo + + !Beta + ispin1=2 + ispin2=2 + + call diexcitation(i,j,k,l,ispin1,ispin2,ref_bitmask,key_out,i_ok,N_int) + + e_double =diag_H_mat_elem(key_out,N_int) + delta_e = 1.d0 / (ref_energy - e_double) + + i_count_connected = 0 + + do a=2, ndet + call get_excitation_degree(key_out,psi_CIS(1,1,a),degree,N_int) + + if (degree /= 1 .and.degree /= 2) cycle + call i_H_j(key_out,psi_CIS(1,1,a),N_int,hij) + if (dabs(hij).le.1.d-10) cycle + + i_count_connected = i_count_connected +1 + + hij_array(i_count_connected) = hij + hij_index(i_count_connected) = a + + enddo + + do a=1,i_count_connected + delta_hij = hij_array(a)*hij_array(a)*delta_e + + delta_H_matrix(hij_index(a),hij_index(a))=delta_H_matrix(hij_index(a),hij_index(a))+delta_hij + + do b=a+1,i_count_connected + delta_hij = hij_array(a)*hij_array(b)*delta_e + + delta_H_matrix(hij_index(a),hij_index(b))=delta_H_matrix(hij_index(a),hij_index(b))+delta_hij + delta_H_matrix(hij_index(b),hij_index(a))=delta_H_matrix(hij_index(b),hij_index(a))+delta_hij + + enddo + enddo + + enddo + enddo + + !different spin contribution + do j=n_core_cis+1,elec_beta_num + + do l=elec_beta_num+1,n_act_cis + ispin1=2 + ispin2=1 + + hij_array=0.d0 + + + call diexcitation(i,j,k,l,ispin1,ispin2,ref_bitmask,key_out,i_ok,N_int) + + e_double =diag_H_mat_elem(key_out,N_int) + delta_e = 1.d0 / (ref_energy - e_double) + + i_count_connected = 0 + + do a=2, ndet + call get_excitation_degree(key_out,psi_CIS(1,1,a),degree,N_int) + + if (degree /= 1 .and.degree /= 2) cycle + + call i_H_j(key_out,psi_CIS(1,1,a),N_int,hij) + + if (dabs(hij).le.1.d-10) cycle + + i_count_connected = i_count_connected +1 + + hij_array(i_count_connected) = hij + hij_index(i_count_connected) = a + + enddo + + do a=1,i_count_connected + delta_hij = hij_array(a)*hij_array(a)*delta_e + + delta_H_matrix(hij_index(a),hij_index(a))=delta_H_matrix(hij_index(a),hij_index(a))+delta_hij + + do b=a+1,i_count_connected + delta_hij = hij_array(a)*hij_array(b)*delta_e + + delta_H_matrix(hij_index(a),hij_index(b))=delta_H_matrix(hij_index(a),hij_index(b))+delta_hij + delta_H_matrix(hij_index(b),hij_index(a))=delta_H_matrix(hij_index(b),hij_index(a))+delta_hij + + enddo + enddo + + enddo + enddo + enddo + enddo + + else !2x2 Matrix diagonalisation + + + delta_H_matrix = 0.d0 + + do i=n_core_cis+1,elec_alpha_num + + do k=elec_alpha_num+1,n_act_cis + + !same spin contribution + do j=i+1,elec_alpha_num + + do l=k+1,n_act_cis + !Alpha + ispin1=1 + ispin2=1 + + call diexcitation(i,j,k,l,ispin1,ispin2,ref_bitmask,key_out,i_ok,N_int) + + e_double =diag_H_mat_elem(key_out,N_int) + + i_count_connected = 0 + accu_H_md=0.d0 + + do a=2, ndet + call get_excitation_degree(key_out,psi_CIS(1,1,a),degree,N_int) + + if (degree /= 1 .and.degree /= 2) cycle + + call i_H_j(key_out,psi_CIS(1,1,a),N_int,hij) + + if (dabs(hij).le.1.d-10) cycle + + i_count_connected = i_count_connected +1 + + hij_array(i_count_connected) = hij + hij_index(i_count_connected) = a + + accu_H_mD=accu_H_mD+(coefs(a))*hij + + enddo + + hmD=accu_H_mD + +! if (dabs(hmD).le.1.d-7) cycle + + epsilon_D=0.5*((e_double-ref_energy)-sqrt((ref_energy -e_double)*(ref_energy -e_double)+4*hmD*hmD)) + if (dabs(hmD*hmD).le.1.d-10.or.dabs(epsilon_D).le.1.d-10) cycle + accu += epsilon_D + accu_2 += hmD*hmD + + delta_e = epsilon_D/(hmD*hmD) + do a=1,i_count_connected + delta_hij = hij_array(a)*hij_array(a) * delta_e + + delta_H_matrix(hij_index(a),hij_index(a))=delta_H_matrix(hij_index(a),hij_index(a))+delta_hij + + do b=a+1,i_count_connected + delta_hij = hij_array(a)*hij_array(b) * delta_e + + delta_H_matrix(hij_index(a),hij_index(b))=delta_H_matrix(hij_index(a),hij_index(b))+delta_hij + delta_H_matrix(hij_index(b),hij_index(a))=delta_H_matrix(hij_index(b),hij_index(a))+delta_hij + + enddo + enddo + + !Beta + ispin1=2 + ispin2=2 + + call diexcitation(i,j,k,l,ispin1,ispin2,ref_bitmask,key_out,i_ok,N_int) + + e_double =diag_H_mat_elem(key_out,N_int) + + i_count_connected = 0 + accu_H_md=0.d0 + + do a=2, ndet + call get_excitation_degree(key_out,psi_CIS(1,1,a),degree,N_int) + + if (degree /= 1 .and.degree /= 2) cycle + + call i_H_j(key_out,psi_CIS(1,1,a),N_int,hij) + + if (dabs(hij).le.1.d-10) cycle + + i_count_connected = i_count_connected +1 + + hij_array(i_count_connected) = hij + hij_index(i_count_connected) = a + + accu_H_mD=accu_H_mD+(coefs(a))*hij + + enddo + + hmD=accu_H_mD + +! if (dabs(hmD).le.1.d-7) cycle + + epsilon_D=0.5*((e_double-ref_energy)-sqrt((ref_energy -e_double)*(ref_energy -e_double)+4*hmD*hmD)) + if (dabs(hmD*hmD).le.1.d-10.or.dabs(epsilon_D).le.1.d-10) cycle + accu += epsilon_D + accu_2 += hmD*hmD + + delta_e = epsilon_D/(hmD*hmD) + + do a=1,i_count_connected + delta_hij = hij_array(a)*hij_array(a) * delta_e + + delta_H_matrix(hij_index(a),hij_index(a))=delta_H_matrix(hij_index(a),hij_index(a))+delta_hij + + do b=a+1,i_count_connected + delta_hij = hij_array(a)*hij_array(b) * delta_e + + delta_H_matrix(hij_index(a),hij_index(b))=delta_H_matrix(hij_index(a),hij_index(b))+delta_hij + delta_H_matrix(hij_index(b),hij_index(a))=delta_H_matrix(hij_index(b),hij_index(a))+delta_hij + + enddo + enddo + + enddo + enddo + + !different spin contribution + do j=n_core_cis+1,elec_beta_num + + do l=elec_beta_num+1,n_act_cis + ispin1=2 + ispin2=1 + + call diexcitation(i,j,k,l,ispin1,ispin2,ref_bitmask,key_out,i_ok,N_int) + + e_double =diag_H_mat_elem(key_out,N_int) + + i_count_connected = 0 + accu_H_md=0.d0 + + do a=2, ndet + call get_excitation_degree(key_out,psi_CIS(1,1,a),degree,N_int) + + if (degree /= 1 .and.degree /= 2) cycle + + call i_H_j(key_out,psi_CIS(1,1,a),N_int,hij) + + if (dabs(hij).le.1.d-10) cycle + + i_count_connected = i_count_connected +1 + + hij_array(i_count_connected) = hij + hij_index(i_count_connected) = a + + accu_H_mD=accu_H_mD+(coefs(a))*hij + + enddo + + hmD=accu_H_mD + + + epsilon_D=0.5*((e_double-ref_energy)-sqrt((ref_energy -e_double)*(ref_energy -e_double)+4*hmD*hmD)) + if (dabs(hmD*hmD).le.1.d-10.or.dabs(epsilon_D).le.1.d-10) cycle + accu += epsilon_D + accu_2 += hmD*hmD + + delta_e = epsilon_D/(hmD*hmD) +! e_double =1.d0/(ref_energy - H_jj_total_diff(key_out,ref_bitmask,HF_energy,N_int)) +! if(dabs(e_double - delta_e)/(dabs(e_double)).gt.1.d-3)then +! print*,delta_e,e_double +! endif + + do a=1,i_count_connected + delta_hij = hij_array(a)*hij_array(a) * delta_e + + delta_H_matrix(hij_index(a),hij_index(a))=delta_H_matrix(hij_index(a),hij_index(a))+delta_hij + + do b=a+1,i_count_connected + delta_hij = hij_array(a)*hij_array(b) * delta_e + + delta_H_matrix(hij_index(a),hij_index(b))=delta_H_matrix(hij_index(a),hij_index(b))+delta_hij + delta_H_matrix(hij_index(b),hij_index(a))=delta_H_matrix(hij_index(b),hij_index(a))+delta_hij + + enddo + enddo + + enddo + enddo + enddo + enddo + + + end if + + end + + subroutine dress_T_con(ref_energy,delta_H_trip,ndet) + + BEGIN_DOC + !Generating all the Triples and, in the loops, the connected Singles + END_DOC + + implicit none + double precision :: e_i,e_r,e_j,e_s,e_k,e_t !epsilons of occupied and virtuals + double precision :: delta_e_ir,delta_e_irj,delta_e_irjs,delta_e_irjsk + double precision :: delta_e,energy !delta epsilons + double precision :: phi_tmp,phi_doub + double precision :: direct,exchg,get_mo_bielec_integral + double precision :: phi_aaa(n_core_cis+1:elec_alpha_num,elec_alpha_num+1:n_act_cis),phi_aab_alpha(n_core_cis+1:elec_alpha_num,elec_alpha_num+1:n_act_cis),phi_aab_beta(n_core_cis+1:elec_alpha_num,elec_alpha_num+1:n_act_cis),phi_aab(n_core_cis+1:elec_alpha_num,elec_alpha_num+1:n_act_cis) + integer :: i,j,k,r,s,t !holes and particles of the Triples + integer :: occ,vir,occ_tmp,vir_tmp !hole and particle of the Singles + integer :: a,b,c,d !variables of the conected Singels + integer :: m,key_ir,key_js,key_is,key_jr,key_kt,key_ks,key_kr,key_it,key_jt,key_tmp + + integer, intent(in) :: ndet + double precision, intent(in) :: ref_energy + double precision, intent(out) :: delta_H_trip(ndet,ndet) + + delta_H_trip=0.d0 +!do i=5,6 + do i=n_core_cis+1,elec_alpha_num + e_i=diagonal_Fock_matrix_mo(i) + !r=7 + do r=elec_alpha_num+1,n_act_cis + e_r=diagonal_Fock_matrix_mo(r) + delta_e_ir=e_i-e_r + key_ir=psi_CIS_adress(i,r) + do j=i+1,elec_alpha_num + e_j=diagonal_Fock_matrix_mo(j) + delta_e_irj=delta_e_ir+e_j + do s=r+1,n_act_cis + e_s=diagonal_Fock_matrix_mo(s) + delta_e_irjs=delta_e_irj-e_s + !alpha-alpha-alpha + do k=j+1,elec_alpha_num + e_k=diagonal_Fock_matrix_mo(k) + delta_e_irjsk=delta_e_irjs+e_k + do t=s+1,n_act_cis + e_t=diagonal_Fock_matrix_mo(t) + delta_e=delta_e_irjsk-e_t + energy=1.d0/(ref_energy-delta_e) + occ=i + a=j + b=k + vir=r + c=s + d=t + direct=get_mo_bielec_integral(a,b,c,d,mo_integrals_map) + exchg =get_mo_bielec_integral(a,b,d,c,mo_integrals_map) + phi_tmp=direct-exchg + do occ=i,k + if (occ==i) then + a=j + b=k + else if (occ==j) then + a=i + b=k + !else if (occ==k) then + ! a=i + ! b=j + !else + ! cycle + end if + do vir=r,t + if (vir==r) then + c=s + d=t + else if (vir==s) then + c=r + d=t + ! else if (vir==t) then + ! c=r + ! d=t + ! else + ! cycle + end if + key_tmp=psi_CIS_adress(occ,vir) + direct=get_mo_bielec_integral(a,b,c,d,mo_integrals_map) + exchg =get_mo_bielec_integral(a,b,d,c,mo_integrals_map) + if (occ==j .and. vir==s .or. occ==k .and. vir==t ) then + phi_aaa(occ,vir)=direct-exchg + else if (occ==i .and. vir==r) then + phi_aaa(occ,vir)=0.d0 + !else if (occ/=i .or. occ/=j .or. occ/=k .or. vir/=r .or. vir/=s .or. vir /=t) then + ! phi_aaa(occ,vir)=0.d0 + else + phi_aaa(occ,vir)=-direct+exchg + endif + phi_doub=phi_tmp*phi_aaa(occ,vir)*energy + delta_H_trip(key_ir,key_tmp)=delta_H_trip(key_ir,key_tmp)+phi_doub + delta_H_trip(key_ir+1,key_tmp+1)=delta_H_trip(key_ir+1,key_tmp+1)+phi_doub + + ! if (phi_doub.ge.1.d-6) then + ! print*,'alpha alpha' + ! print*,'occ,vir,key_tmp',occ,vir,key_tmp + ! endif + + ! if (key_ir==key_tmp) then + ! print*,'ir=tmp' + ! print*,'i,r,j,s,k,t',i,r,j,s,k,t + ! print*,'occ,vir',occ,vir + ! !stop + ! end if + + enddo + enddo + enddo + enddo + !alpha-alpha-beta + do k=n_core_cis+1,elec_alpha_num + e_k=diagonal_Fock_matrix_mo(k) + delta_e_irjsk=delta_e_irjs+e_k + do t=elec_alpha_num+1,n_act_cis + e_t=diagonal_Fock_matrix_mo(t) + delta_e=delta_e_irjsk-e_t + energy=1.d0/(ref_energy-delta_e) + occ=i + a=j + b=k + vir=r + c=s + d=t + direct=get_mo_bielec_integral(a,b,c,d,mo_integrals_map) + phi_tmp=direct-exchg + do occ=i,k + if (occ==i) then + a=j + b=k + else if (occ==j) then + a=i + b=k + else if (occ==k) then + a=i + b=j + else + cycle + end if + do vir=r,t + if (vir==r) then + c=s + d=t + else if (vir==s) then + c=r + d=t + else if (vir==t) then + c=r + d=s + else + cycle + end if + key_tmp=psi_CIS_adress(occ,vir) + direct=get_mo_bielec_integral(a,b,c,d,mo_integrals_map) + if (occ==j .and. vir==s .or. occ==k .and. vir==t ) then + phi_aab(occ,vir)=direct + else if (occ==i .and. vir==r .or. occ==i .and. vir==t .or. occ==j .and. vir==t .or. occ==k .and. vir==r .or. occ==k .and. vir==s) then + phi_aab(occ,vir)=0.d0 + !else if (occ/=i .or. occ/=j .or. occ/=k .or. vir/=r .or. vir/=s .or. vir /=t) then + ! phi_aab(occ,vir)=0.d0 + + else + phi_aab(occ,vir)=-direct + endif + + phi_doub=phi_tmp*phi_aab(occ,vir)*energy + delta_H_trip(key_ir,key_tmp+1)=delta_H_trip(key_ir,key_tmp+1)+phi_doub + delta_H_trip(key_ir+1,key_tmp)=delta_H_trip(key_ir,key_tmp+1)+phi_doub + ! if (phi_doub.ge.1.d-6) then + ! if (key_ir==key_tmp) then + ! print*,'ir=tmp beta ' + ! print*,'i,r,j,s,k,t',i,r,j,s,k,t + ! print*,'occ,vir',occ,vir + ! !stop + ! end if + ! end if + + ! if (phi_doub.ge.1.d-6) then + ! if (vir==r .or. vir==s) then + ! print*,'alpha' + ! print*,'occ,vir,key_tmp',occ,vir,key_tmp + + ! else if (vir==t) then + ! print*,'beta' + ! print*,'occ,vir,key_tmp',occ,vir,key_tmp + + ! endif + ! print*,'vir=?',vir + ! endif + + enddo + enddo + enddo + enddo + enddo + enddo + enddo + enddo + +!delta_H_trip=0.d0 + + + +!!generating the Singles included in the Triples +!! do m=1,n_state_CIS +!do i=n_core_cis+1,elec_alpha_num +! e_i=diagonal_Fock_matrix_mo(i) + + + +!!do r=elec_alpha_num+1,n_state_CIS +! do r=elec_alpha_num+1,n_act_cis +! e_r=diagonal_Fock_matrix_mo(r) +! delta_e_ir=e_i-e_r +!key_ir=psi_CIS_adress(i,r) + +! !alpha-alpha-x (=beta-beta-x) +! do j=i+1,elec_alpha_num +! e_j=diagonal_Fock_matrix_mo(j) +! delta_e_irj=delta_e_ir+e_j +!key_jr=psi_CIS_adress(j,r) +! do s=r+1,n_act_cis +! e_s=diagonal_Fock_matrix_mo(s) +! delta_e_irjs=delta_e_irj-e_s + +! key_is=psi_CIS_adress(i,s) +! key_js=psi_CIS_adress(j,s) + + +! !alpha-alpha-alpha (=beta-beta-beta) +! do k=j+1,elec_alpha_num +! e_k=diagonal_Fock_matrix_mo(k) +! delta_e_irjsk=delta_e_irjs+e_k +! key_kr=psi_CIS_adress(k,r) +! key_ks=psi_CIS_adress(k,s) + +! do t=s+1,n_act_cis +! e_t=diagonal_Fock_matrix_mo(t) +! delta_e=delta_e_irjsk-e_t + +! key_it=psi_CIS_adress(i,t) +! key_jt=psi_CIS_adress(j,t) +! key_kt=psi_CIS_adress(k,t) + + + +! !generating the connected singles +! do occ=i,k +! if (occ==i) then +! a=j +! b=k +! else if (occ==j) then +! a=i +! b=k +! else if (occ==k) then +! a=i +! b=j +! else +! cycle +! end if + +! do vir=r,t +! if (vir==r) then +! c=s +! d=t +! else if (vir==s) then +! c=r +! d=t +! else if (vir==t) then +! c=r +! d=s +! else +! cycle +! end if + +! if (occ==i .and. vir==r .or. occ==j .and. vir==s .or. occ==k .and. vir==t ) then +! direct=get_mo_bielec_integral(a,b,c,d,mo_integrals_map) +! exchg =get_mo_bielec_integral(a,b,d,c,mo_integrals_map) + +! phi_aaa(occ,vir)=direct-exchg +! else +! direct=get_mo_bielec_integral(a,b,c,d,mo_integrals_map) +! exchg =get_mo_bielec_integral(a,b,d,c,mo_integrals_map) + +! phi_aaa(occ,vir)=-direct+exchg +! endif + +! ! phi_aaa(occ,vir)=phi_aaa(occ,vir)+phi_tmp + +! enddo +! enddo + +! phi_tmp=(phi_aaa(i,r)*phi_aaa(j,s))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_ir,key_js)=delta_H_trip(key_ir,key_js)+phi_tmp +! delta_H_trip(key_ir+1,key_js+1)=delta_H_trip(key_ir+1,key_js+1)+phi_tmp + +! phi_tmp=(phi_aaa(i,r)*phi_aaa(k,t))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_ir,key_kt)=delta_H_trip(key_ir,key_kt)+phi_tmp +! delta_H_trip(key_ir+1,key_kt+1)=delta_H_trip(key_ir+1,key_kt+1)+phi_tmp + +! phi_tmp=(phi_aaa(i,r)*phi_aaa(j,t))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_ir,key_jt)=delta_H_trip(key_ir,key_jt)+phi_tmp +! delta_H_trip(key_ir+1,key_jt+1)=delta_H_trip(key_ir+1,key_jt+1)+phi_tmp + +! phi_tmp=(phi_aaa(i,r)*phi_aaa(k,s))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_ir,key_ks)=delta_H_trip(key_ir,key_ks)+phi_tmp +! delta_H_trip(key_ir+1,key_ks+1)=delta_H_trip(key_ir+1,key_ks+1)+phi_tmp + + +! phi_tmp=(phi_aaa(i,s)*phi_aaa(j,r))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_is,key_jr)=delta_H_trip(key_is,key_jr)+phi_tmp +! delta_H_trip(key_is+1,key_jr+1)=delta_H_trip(key_is+1,key_jr+1)+phi_tmp + +! phi_tmp=(phi_aaa(i,s)*phi_aaa(k,t))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_is,key_kt)=delta_H_trip(key_is,key_kt)+phi_tmp +! delta_H_trip(key_is+1,key_kt+1)=delta_H_trip(key_is+1,key_kt+1)+phi_tmp + +! phi_tmp=(phi_aaa(i,s)*phi_aaa(j,t))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_is,key_jt)=delta_H_trip(key_is,key_jt)+phi_tmp +! delta_H_trip(key_is+1,key_jt+1)=delta_H_trip(key_is+1,key_jt+1)+phi_tmp + +! phi_tmp=(phi_aaa(i,s)*phi_aaa(k,r))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_is,key_kr)=delta_H_trip(key_is,key_kr)+phi_tmp +! delta_H_trip(key_is+1,key_kr+1)=delta_H_trip(key_is+1,key_kr+1)+phi_tmp + +! phi_tmp=(phi_aaa(i,t)*phi_aaa(j,s))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_it,key_js)=delta_H_trip(key_it,key_js)+phi_tmp +! delta_H_trip(key_it+1,key_js+1)=delta_H_trip(key_it+1,key_js+1)+phi_tmp + +! phi_tmp=(phi_aaa(i,t)*phi_aaa(j,r))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_it,key_jr)=delta_H_trip(key_it,key_jr)+phi_tmp +! delta_H_trip(key_it+1,key_jr+1)=delta_H_trip(key_it+1,key_jr+1)+phi_tmp + +! phi_tmp=(phi_aaa(i,t)*phi_aaa(k,s))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_it,key_ks)=delta_H_trip(key_it,key_ks)+phi_tmp +! delta_H_trip(key_it+1,key_ks+1)=delta_H_trip(key_it+1,key_ks+1)+phi_tmp + +! phi_tmp=(phi_aaa(i,t)*phi_aaa(k,r))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_it,key_kr)=delta_H_trip(key_it,key_kr)+phi_tmp +! delta_H_trip(key_it+1,key_kr+1)=delta_H_trip(key_it+1,key_kr+1)+phi_tmp + + +! phi_tmp=(phi_aaa(j,r)*phi_aaa(k,t))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_jr,key_kt)=delta_H_trip(key_jr,key_kt)+phi_tmp +! delta_H_trip(key_jr+1,key_kt+1)=delta_H_trip(key_jr+1,key_kt+1)+phi_tmp + +! phi_tmp=(phi_aaa(j,r)*phi_aaa(k,s))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_jr,key_ks)=delta_H_trip(key_jr,key_ks)+phi_tmp +! delta_H_trip(key_jr+1,key_ks+1)=delta_H_trip(key_jr+1,key_ks+1)+phi_tmp + +! phi_tmp=(phi_aaa(j,s)*phi_aaa(k,t))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_js,key_kt)=delta_H_trip(key_js,key_kt)+phi_tmp +! delta_H_trip(key_js+1,key_kt+1)=delta_H_trip(key_js+1,key_kt+1)+phi_tmp + +! phi_tmp=(phi_aaa(j,s)*phi_aaa(k,r))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_js,key_kr)=delta_H_trip(key_js,key_kr)+phi_tmp +! delta_H_trip(key_js+1,key_kr+1)=delta_H_trip(key_js+1,key_kr+1)+phi_tmp + +! phi_tmp=(phi_aaa(j,t)*phi_aaa(k,s))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_jt,key_ks)=delta_H_trip(key_jt,key_ks)+phi_tmp +! delta_H_trip(key_jt+1,key_ks+1)=delta_H_trip(key_jt+1,key_ks+1)+phi_tmp + +! phi_tmp=(phi_aaa(j,t)*phi_aaa(k,r))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_jt,key_kr)=delta_H_trip(key_jt,key_kr)+phi_tmp +! delta_H_trip(key_jt+1,key_kr+1)=delta_H_trip(key_jt+1,key_kr+1)+phi_tmp + +! enddo +! enddo + +! !alpha-alpha-beta (=beta-beta-alpha) +! do k=n_core_cis+1,elec_beta_num +! e_k=diagonal_Fock_matrix_mo(k) +! delta_e_irjsk=delta_e_irjs+e_k + +! key_kr=psi_CIS_adress(k,r) +! key_ks=psi_CIS_adress(k,s) + + + +! do t=elec_beta_num+1,n_act_cis +! e_t=diagonal_Fock_matrix_mo(t) +! delta_e=delta_e_irjsk-e_t + +! key_it=psi_CIS_adress(i,t) +! key_jt=psi_CIS_adress(j,t) +! key_kt=psi_CIS_adress(k,t) + + +! !generating the connected singles alpha (beta) +! do occ=i,j +! if (occ==i) then +! a=j +! b=k +! else if (occ==j) then +! a=i +! b=k +! end if + +! do vir = r,s +! if (vir==r) then +! c=s +! d=t +! else if (vir==s) then +! c=r +! d=t +! end if + +! direct=get_mo_bielec_integral(a,b,c,d,mo_integrals_map) + +! if (occ==i .and. vir==r .or. occ==j .and. vir==s) then +! phi_aab_alpha(occ,vir)=direct +! else if (occ==i .and. vir==s .or. occ==j .and. vir==r) then +! phi_aab_alpha(occ,vir)=-1.d0*direct + +! end if + +!! phi_aab_alpha(occ,vir)=phi_aab_alpha(occ,vir)+phi_tmp + +! enddo +! enddo + +! !alpha single alpha single connected +! phi_tmp=(phi_aab_alpha(i,r)*phi_aab_alpha(j,s))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_ir,key_js)=delta_H_trip(key_ir,key_js)+phi_tmp +! delta_H_trip(key_ir+1,key_js+1)=delta_H_trip(key_ir+1,key_js+1)+phi_tmp + + +! ! phi_tmp=(phi_aab_alpha(i,r)*phi_aab_alpha(i,s))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! ! delta_H_trip(key_ir,key_is)=delta_H_trip(key_ir,key_is)+phi_tmp +! ! delta_H_trip(key_ir+1,key_is+1)=delta_H_trip(key_ir+1,key_is+1)+phi_tmp + + +! phi_tmp=(phi_aab_alpha(i,s)*phi_aab_alpha(j,r))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_is,key_jr)=delta_H_trip(key_is,key_jr)+phi_tmp +! delta_H_trip(key_is+1,key_jr+1)=delta_H_trip(key_is+1,key_jr+1)+phi_tmp + + + + + +! !generating the conncected singles beta (alpha) +! ! do occ=i,k +! ! if (occ==i) then +! ! a=j +! ! b=k +! ! else if (occ==j) then +! ! a=i +! ! b=k +! ! else if (occ==k) then +! ! +! ! a=i +! ! b=j +! ! end if +! ! do vir=r,t +! ! if (vir==r) then +! ! c=s +! ! d=t +! ! else if (vir==s) then +! ! c=r +! ! d=t +! ! else if (vir==k) then +! ! +! ! c=r +! ! d=t +! ! end if + +! occ=k +! a=i +! b=j +! vir=t +! c=r +! d=s +! direct=get_mo_bielec_integral(a,b,c,d,mo_integrals_map) +! exchg =get_mo_bielec_integral(a,b,d,c,mo_integrals_map) + +! phi_aab_beta(occ,vir)=direct-exchg +! ! enddo +! ! enddo +!! phi_aab_beta(occ,vir)=phi_aab_beta(occ,vir)+phi_tmp +! + + + + +! !alpha single beta single connected +! phi_tmp=(phi_aab_alpha(i,r)*phi_aab_beta(k,t))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_ir,key_kt+1)=delta_H_trip(key_ir,key_kt+1)+phi_tmp +! delta_H_trip(key_ir+1,key_kt)=delta_H_trip(key_ir+1,key_kt)+phi_tmp + +! phi_tmp=(phi_aab_alpha(i,s)*phi_aab_beta(k,t))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_is,key_kt+1)=delta_H_trip(key_is,key_kt+1)+phi_tmp +! delta_H_trip(key_is+1,key_kt)=delta_H_trip(key_is+1,key_kt)+phi_tmp + + +! phi_tmp=(phi_aab_alpha(j,r)*phi_aab_beta(k,t))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_jr,key_kt+1)=delta_H_trip(key_jr,key_kt+1)+phi_tmp +! delta_H_trip(key_jr+1,key_kt)=delta_H_trip(key_jr+1,key_kt)+phi_tmp + + +! phi_tmp=(phi_aab_alpha(j,s)*phi_aab_beta(k,t))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_js+1,key_kt)=delta_H_trip(key_js+1,key_kt)+phi_tmp +! delta_H_trip(key_js,key_kt+1)=delta_H_trip(key_js,key_kt+1)+phi_tmp + + +! enddo +! enddo +! enddo +! enddo +! enddo +!enddo + + END + +subroutine diexcitation(i,j,k,l,ispin1,ispin2,key_in,key_out,i_ok,Nint) + implicit none + use bitmasks + ! realize the double excitation i-->k (ispin1) + j-->l (ispin2) on key_in + ! returns key_out and i_ok (i_ok = 0 means not possible, i_ok = 1 means the excitation was possible) + integer, intent(in) :: ispin1,ispin2,i,j,k,l,Nint + integer(bit_kind), intent(in) :: key_in(Nint,2) + integer, intent(out):: i_ok + integer(bit_kind), intent(out):: key_out(Nint,2) + integer :: k_hole,j_hole,k_particl,j_particl,i_nint,Nelec_alpha,Nelec_beta + integer :: i_test_hole,i_test_particl + key_out = key_in + + k_hole = ishft(i-1,-5)+1 + j_hole = i-ishft(k_hole-1,5)-1 + i_test_hole = ibset(0,j_hole) + if(iand(key_in(k_hole,ispin1),i_test_hole).ne.i_test_hole)then + i_ok = 0 + return + endif + key_out(k_hole,ispin1) = ibclr(key_out(k_hole,ispin1),j_hole) + k_particl = ishft(k-1,-5)+1 + j_particl = k-ishft(k_particl-1,5)-1 + i_test_particl= ibset(0,j_particl) + if(iand(key_in(k_particl,ispin1),i_test_particl).ne.0)then + i_ok = 0 + return + endif + key_out(k_particl,ispin1) = ibset(key_out(k_particl,ispin1),j_particl) + + k_hole = ishft(j-1,-5)+1 + j_hole = j-ishft(k_hole-1,5)-1 + i_test_hole = ibset(0,j_hole) + if(iand(key_in(k_hole,ispin2),i_test_hole).ne.i_test_hole)then + i_ok = 0 + return + endif + key_out(k_hole,ispin2) = ibclr(key_out(k_hole,ispin2),j_hole) + k_particl = ishft(l-1,-5)+1 + j_particl = l-ishft(k_particl-1,5)-1 + i_test_particl = ibset(0,j_particl) + if(iand(key_in(k_particl,ispin2),i_test_particl).ne.0)then + i_ok = 0 + return + endif + key_out(k_particl,ispin2) = ibset(key_out(k_particl,ispin2),j_particl) + i_ok = 1 + end + diff --git a/src/CIS_dressed/README.rst b/src/CIS_dressed/README.rst new file mode 100644 index 00000000..bb6eeeea --- /dev/null +++ b/src/CIS_dressed/README.rst @@ -0,0 +1,28 @@ +Needed Modules +============== + +.. Do not edit this section. It was auto-generated from the +.. NEEDED_MODULES file. + +* `AOs `_ +* `BiInts `_ +* `Bitmask `_ +* `Dets `_ +* `Electrons `_ +* `Ezfio_files `_ +* `Hartree_Fock `_ +* `MonoInts `_ +* `MOs `_ +* `Nuclei `_ +* `Output `_ +* `Utils `_ +* `DensityMatrix `_ + +Documentation +============= + +.. Do not edit this section. It was auto-generated from the +.. NEEDED_MODULES file. + + + diff --git a/src/CIS_dressed/cis_dressed.ezfio_config b/src/CIS_dressed/cis_dressed.ezfio_config new file mode 100644 index 00000000..fd9f291a --- /dev/null +++ b/src/CIS_dressed/cis_dressed.ezfio_config @@ -0,0 +1,6 @@ +cis_dressed + n_state_cis integer + n_core_cis integer + n_act_cis integer + mp2_dressing logical + standard_doubles logical diff --git a/src/CIS_dressed/density_matrix_suroutine.irp.f b/src/CIS_dressed/density_matrix_suroutine.irp.f new file mode 100644 index 00000000..59e099a8 --- /dev/null +++ b/src/CIS_dressed/density_matrix_suroutine.irp.f @@ -0,0 +1,66 @@ + + subroutine get_dm_from_psi(dets_in,u_in,sze,dim_in,Nint,dm_alpha,dm_beta) + implicit none + BEGIN_DOC + ! Alpha and beta one-body density matrix + ! + ! dets_in :: bitsrings corresponding to the determinants in the wave function + ! + ! u_in :: coefficients of the wave function + ! + ! sze :: number of determinants in the wave function + ! + ! dim_in :: physical dimension of the array u_in and dets_in + ! + ! Nint :: should be equal to N_int + ! + ! dm_alpha :: alpha one body density matrix + ! + ! dm_beta :: beta one body density matrix + END_DOC + use bitmasks + integer, intent(in) :: sze,dim_in,Nint + integer(bit_kind), intent(in) :: dets_in(Nint,2,dim_in) + double precision, intent(in) :: u_in(dim_in) + double precision, intent(out) :: dm_alpha(mo_tot_num,mo_tot_num) + double precision, intent(out) :: dm_beta(mo_tot_num,mo_tot_num) + + integer :: j,k,l + integer :: occ(N_int*bit_kind_size,2) + double precision :: ck, cl, ckl + double precision :: phase + integer :: h1,h2,p1,p2,s1,s2, degree + integer :: exc(0:2,2,2),n_occ_alpha + dm_alpha = 0.d0 + dm_beta = 0.d0 + + do k=1,sze + call bitstring_to_list(dets_in(1,1,k), occ(1,1), n_occ_alpha, N_int) + call bitstring_to_list(dets_in(1,2,k), occ(1,2), n_occ_alpha, N_int) + ck = u_in(k) + do l=1,elec_alpha_num + j = occ(l,1) + dm_alpha(j,j) += ck*ck + enddo + do l=1,elec_beta_num + j = occ(l,2) + dm_beta(j,j) += ck*ck + enddo + do l=1,k-1 + call get_excitation_degree(dets_in(1,1,k),dets_in(1,1,l),degree,N_int) + if (degree /= 1) then + cycle + endif + call get_mono_excitation(dets_in(1,1,k),dets_in(1,1,l),exc,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + ckl = ck * u_in(l) * phase + if (s1==1) then + dm_alpha(h1,p1) += ckl + dm_alpha(p1,h1) += ckl + else + dm_beta(h1,p1) += ckl + dm_beta(p1,h1) += ckl + endif + enddo + enddo + end diff --git a/src/CIS_dressed/natural_particl_orbitals.irp.f b/src/CIS_dressed/natural_particl_orbitals.irp.f new file mode 100644 index 00000000..532f640f --- /dev/null +++ b/src/CIS_dressed/natural_particl_orbitals.irp.f @@ -0,0 +1,92 @@ + BEGIN_PROVIDER [double precision, particle_natural_orb_CIS_properties, (6,n_state_cis)] +&BEGIN_PROVIDER [double precision, CIS_states_properties, (6,n_state_cis)] + implicit none + BEGIN_DOC +! properties of the natural orbital of the particle of the various n_state_cis eigenvectors of the CIS matrix +! +! You first build the density matrix of the one eigenvector and you take off the Hartree Fock density matrix +! +! particl(i,j)(state = k) == dm(i,j)(Hartree Fock) - dm(i,j)(state = k) +! +! you diagonalize particl(i,j) and the first eigenvector is the natural orbital corresponding to the particl +! +! that is specific to the excitation in the CIS state +! +! particle_natural_orb_CIS_properties(i,1) = +! +! particle_natural_orb_CIS_properties(i,2) = +! +! particle_natural_orb_CIS_properties(i,3) = +! +! particle_natural_orb_CIS_properties(i,5) = +! +! particle_natural_orb_CIS_properties(i,6) = +! +! particle_natural_orb_CIS_properties(i,7) = +! +! CIS_states_properties(i,1:6) = the same but for the hole state i + END_DOC + + integer :: i,j,k,l + double precision :: dm_alpha(mo_tot_num,mo_tot_num) + double precision :: dm_beta(mo_tot_num,mo_tot_num) + double precision :: dm(mo_tot_num,mo_tot_num) + double precision :: eigvalues(mo_tot_num) + double precision :: eigvectors(mo_tot_num,mo_tot_num) + double precision :: accu_n_elec,c_k + + + do i = 1, n_state_cis + print*,' state cis = ',i + call get_dm_from_psi(psi_CIS,coefs_CIS(1,i),size_psi_CIS,size_psi_CIS,N_int,dm_alpha,dm_beta) + + dm = dm_alpha + dm_beta + call get_properties_from_density_matrix(dm,CIS_states_properties(1,i)) + dm = -dm + do k = 1, elec_alpha_num + dm(k,k) += 1.d0 + enddo + do k = 1, elec_beta_num + dm(k,k) += 1.d0 + enddo + call lapack_diag(eigvalues,eigvectors,dm,mo_tot_num,mo_tot_num) + accu_n_elec = 0.d0 + do k = 1, mo_tot_num + accu_n_elec += eigvalues(k) + enddo + do k = 1, mo_tot_num + do l = 1, mo_tot_num + c_k = eigvectors(k,j) * eigvectors(l,j) + particle_natural_orb_CIS_properties(1,i) += c_k * mo_dipole_x(k,l) + particle_natural_orb_CIS_properties(2,i) += c_k * mo_dipole_y(k,l) + particle_natural_orb_CIS_properties(3,i) += c_k * mo_dipole_z(k,l) + particle_natural_orb_CIS_properties(4,i) += c_k * mo_spread_x(k,l) + particle_natural_orb_CIS_properties(5,i) += c_k * mo_spread_y(k,l) + particle_natural_orb_CIS_properties(6,i) += c_k * mo_spread_z(k,l) + enddo + enddo + enddo + +END_PROVIDER +subroutine get_properties_from_density_matrix(dm,properties) + implicit none + double precision, intent(in) :: dm(mo_tot_num,mo_tot_num) + double precision, intent(out) :: properties(6) + integer :: k,l + double precision :: c_k + do k = 1, 6 + properties(k) = 0.d0 + enddo + do k = 1, mo_tot_num + do l = 1, mo_tot_num + c_k = dm(k,l) + properties(1) += c_k * mo_dipole_x(k,l) + properties(2) += c_k * mo_dipole_y(k,l) + properties(3) += c_k * mo_dipole_z(k,l) + properties(4) += c_k * mo_spread_x(k,l) + properties(5) += c_k * mo_spread_y(k,l) + properties(6) += c_k * mo_spread_z(k,l) + enddo + enddo + +end diff --git a/src/CIS_dressed/options.irp.f b/src/CIS_dressed/options.irp.f new file mode 100644 index 00000000..59b67c4f --- /dev/null +++ b/src/CIS_dressed/options.irp.f @@ -0,0 +1,90 @@ +BEGIN_PROVIDER [ integer , n_state_cis ] + implicit none + BEGIN_DOC + ! Number of states asked for the CIS vector + END_DOC + + logical :: has + PROVIDE ezfio_filename + call ezfio_has_cis_dressed_n_state_cis(has) + if (has) then + call ezfio_get_cis_dressed_n_state_cis(n_state_cis) + else + n_state_cis = 5 + call ezfio_set_cis_dressed_n_state_cis(n_state_cis) + endif + +END_PROVIDER + +BEGIN_PROVIDER [ integer , n_core_cis] + implicit none + BEGIN_DOC + ! Number of states asked for the CIS vector + END_DOC + + logical :: has + PROVIDE ezfio_filename + call ezfio_has_cis_dressed_n_core_cis(has) + if (has) then + call ezfio_get_cis_dressed_n_core_cis(n_core_cis) + else + n_core_cis = 0 + call ezfio_set_cis_dressed_n_core_cis(n_core_cis) + endif + +END_PROVIDER + +BEGIN_PROVIDER [ integer , n_act_cis] + implicit none + BEGIN_DOC + ! Number of states asked for the CIS vector + END_DOC + + logical :: has + PROVIDE ezfio_filename + call ezfio_has_cis_dressed_n_act_cis(has) + if (has) then + call ezfio_get_cis_dressed_n_act_cis(n_act_cis) + else + n_act_cis = mo_tot_num + call ezfio_set_cis_dressed_n_act_cis(n_act_cis) + endif + +END_PROVIDER + +BEGIN_PROVIDER [ logical , mp2_dressing] + implicit none + BEGIN_DOC + ! Number of states asked for the CIS vector + END_DOC + + logical :: has + PROVIDE ezfio_filename + call ezfio_has_cis_dressed_mp2_dressing(has) + if (has) then + call ezfio_get_cis_dressed_mp2_dressing(mp2_dressing) + else + mp2_dressing = .False. + call ezfio_set_cis_dressed_mp2_dressing(mp2_dressing) + endif + +END_PROVIDER + +BEGIN_PROVIDER [ logical , standard_doubles] + implicit none + BEGIN_DOC + ! Number of states asked for the CIS vector + END_DOC + + logical :: has + PROVIDE ezfio_filename + call ezfio_has_cis_dressed_standard_doubles(has) + if (has) then + call ezfio_get_cis_dressed_standard_doubles(standard_doubles) + else + standard_doubles = .True. + call ezfio_set_cis_dressed_standard_doubles(standard_doubles) + endif + +END_PROVIDER + diff --git a/src/DensityMatrix/README.rst b/src/DensityMatrix/README.rst index d0734ed9..af3c405c 100644 --- a/src/DensityMatrix/README.rst +++ b/src/DensityMatrix/README.rst @@ -8,6 +8,39 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. +`iunit_two_body_dm_aa `_ + Temporary files for 2-body dm calculation + +`iunit_two_body_dm_ab `_ + Temporary files for 2-body dm calculation + +`iunit_two_body_dm_bb `_ + Temporary files for 2-body dm calculation + +`one_body_dm_a `_ + Alpha and beta one-body density matrix + +`one_body_dm_b `_ + Alpha and beta one-body density matrix + +`two_body_dm_diag_aa `_ + diagonal part of the two body density matrix + +`two_body_dm_diag_ab `_ + diagonal part of the two body density matrix + +`two_body_dm_diag_bb `_ + diagonal part of the two body density matrix + +`det_coef_provider `_ + Undocumented + +`det_num `_ + Undocumented + +`det_provider `_ + Undocumented + Needed Modules @@ -28,4 +61,5 @@ Needed Modules * `Nuclei `_ * `Output `_ * `Utils `_ +* `DensityMatrix `_ diff --git a/src/DensityMatrix/det_num.irp.f b/src/DensityMatrix/det_num.irp.f index 7c73b0a4..6ba0c1b4 100644 --- a/src/DensityMatrix/det_num.irp.f +++ b/src/DensityMatrix/det_num.irp.f @@ -43,14 +43,14 @@ -0.316337207748718D-01 & /) - do i=1,10 - call write_bitstring( 6, det_provider(1,1,i), N_int ) - enddo - print *, '' - do i=1,10 - call write_bitstring( 6, det_provider(1,2,i), N_int ) - enddo - print *, '' +!do i=1,10 +! call write_bitstring( 6, det_provider(1,1,i), N_int ) +!enddo +!print *, '' +!do i=1,10 +! call write_bitstring( 6, det_provider(1,2,i), N_int ) +!enddo +!print *, '' END_PROVIDER diff --git a/src/Dets/davidson.irp.f b/src/Dets/davidson.irp.f index f0dc2437..09fad255 100644 --- a/src/Dets/davidson.irp.f +++ b/src/Dets/davidson.irp.f @@ -275,3 +275,4 @@ subroutine davidson_diag(dets_in,u_in,energies,dim_in,sze,N_st,Nint) ) end + diff --git a/src/Dets/s2.irp.f b/src/Dets/s2.irp.f index c970d628..b218e7cb 100644 --- a/src/Dets/s2.irp.f +++ b/src/Dets/s2.irp.f @@ -1,11 +1,12 @@ subroutine get_s2(key_i,key_j,phase,Nint) implicit none + use bitmasks BEGIN_DOC ! Returns END_DOC integer, intent(in) :: Nint - integer, intent(in) :: key_i(Nint,2) - integer, intent(in) :: key_j(Nint,2) + integer(bit_kind), intent(in) :: key_i(Nint,2) + integer(bit_kind), intent(in) :: key_j(Nint,2) double precision, intent(out) :: phase integer :: exc(0:2,2,2) integer :: degree @@ -32,3 +33,38 @@ subroutine get_s2(key_i,key_j,phase,Nint) end select end +BEGIN_PROVIDER [ double precision, S_z ] +&BEGIN_PROVIDER [ double precision, S_z2_Sz ] + implicit none + + S_z = 0.5d0*dble(elec_alpha_num-elec_beta_num) + S_z2_Sz = S_z*(S_z-1.d0) + +END_PROVIDER + + +subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2) + implicit none + use bitmasks + integer(bit_kind), intent(in) :: psi_keys_tmp(N_int,2,nmax) + integer, intent(in) :: n,nmax + double precision, intent(in) :: psi_coefs_tmp(nmax) + double precision, intent(out) :: s2 + integer :: i,j,l + double precision :: s2_tmp + s2 = S_z2_Sz + !$OMP PARALLEL DO DEFAULT(NONE) & + !$OMP PRIVATE(i,j,s2_tmp) SHARED(n,psi_coefs_tmp,psi_keys_tmp,N_int) & + !$OMP REDUCTION(+:s2) SCHEDULE(dynamic) + do i = 1, n + call get_s2(psi_keys_tmp(1,1,i),psi_keys_tmp(1,1,i),s2_tmp,N_int) +! print*,'s2_tmp = ',s2_tmp + do j = 1, n + call get_s2(psi_keys_tmp(1,1,i),psi_keys_tmp(1,1,j),s2_tmp,N_int) + if (s2_tmp == 0.d0) cycle + s2 += psi_coefs_tmp(i)*psi_coefs_tmp(j)*s2_tmp + enddo + enddo + !$OMP END PARALLEL DO +end + diff --git a/src/Hartree_Fock/README.rst b/src/Hartree_Fock/README.rst index bd52e526..f6b194b8 100644 --- a/src/Hartree_Fock/README.rst +++ b/src/Hartree_Fock/README.rst @@ -21,19 +21,19 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`fock_matrix_alpha_ao `_ +`fock_matrix_alpha_ao `_ Alpha Fock matrix in AO basis set -`fock_matrix_alpha_mo `_ +`fock_matrix_alpha_mo `_ Fock matrix on the MO basis -`fock_matrix_beta_ao `_ +`fock_matrix_beta_ao `_ Alpha Fock matrix in AO basis set -`fock_matrix_beta_mo `_ +`fock_matrix_beta_mo `_ Fock matrix on the MO basis -`fock_matrix_diag_mo `_ +`fock_matrix_diag_mo `_ Fock matrix on the MO basis. For open shells, the ROHF Fock Matrix is .br @@ -48,7 +48,7 @@ Documentation K = Fb - Fa .br -`fock_matrix_mo `_ +`fock_matrix_mo `_ Fock matrix on the MO basis. For open shells, the ROHF Fock Matrix is .br @@ -63,49 +63,49 @@ Documentation K = Fb - Fa .br -`hf_energy `_ +`hf_energy `_ Hartree-Fock energy -`hf_density_matrix_ao `_ +`hf_density_matrix_ao `_ Density matrix in the AO basis -`hf_density_matrix_ao_alpha `_ +`hf_density_matrix_ao_alpha `_ Alpha and Beta density matrix in the AO basis -`hf_density_matrix_ao_beta `_ +`hf_density_matrix_ao_beta `_ Alpha and Beta density matrix in the AO basis -`diagonal_fock_matrix_mo `_ +`diagonal_fock_matrix_mo `_ Diagonal Fock matrix in the MO basis -`eigenvectors_fock_matrix_mo `_ +`eigenvectors_fock_matrix_mo `_ Diagonal Fock matrix in the MO basis -`scf_iteration `_ +`scf_iteration `_ Undocumented -`do_diis `_ +`do_diis `_ If True, compute integrals on the fly -`n_it_scf_max `_ +`n_it_scf_max `_ Maximum number of SCF iterations -`thresh_scf `_ +`thresh_scf `_ Threshold on the convergence of the Hartree Fock energy -`bi_elec_ref_bitmask_energy `_ +`bi_elec_ref_bitmask_energy `_ Energy of the reference bitmask used in Slater rules -`kinetic_ref_bitmask_energy `_ +`kinetic_ref_bitmask_energy `_ Energy of the reference bitmask used in Slater rules -`mono_elec_ref_bitmask_energy `_ +`mono_elec_ref_bitmask_energy `_ Energy of the reference bitmask used in Slater rules -`nucl_elec_ref_bitmask_energy `_ +`nucl_elec_ref_bitmask_energy `_ Energy of the reference bitmask used in Slater rules -`ref_bitmask_energy `_ +`ref_bitmask_energy `_ Energy of the reference bitmask used in Slater rules diff --git a/src/Makefile.config.example b/src/Makefile.config.example deleted file mode 100644 index 82e0ca27..00000000 --- a/src/Makefile.config.example +++ /dev/null @@ -1,29 +0,0 @@ -OPENMP =1 -PROFILE =0 -DEBUG = 0 - -IRPF90_FLAGS+= --align=32 -FC = ifort -g -FCFLAGS= -FCFLAGS+= -xHost -#FCFLAGS+= -xAVX -FCFLAGS+= -O2 -FCFLAGS+= -ip -FCFLAGS+= -opt-prefetch -FCFLAGS+= -ftz -MKL=-mkl=parallel - -ifeq ($(PROFILE),1) -FC += -p -g -CXX += -pg -endif - -ifeq ($(OPENMP),1) -FC += -openmp -CXX += -fopenmp -endif - -ifeq ($(DEBUG),1) -FC += -C -traceback -fpe0 -#FCFLAGS =-O0 -endif diff --git a/src/MonoInts/README.rst b/src/MonoInts/README.rst index 03ec9939..e36763d9 100644 --- a/src/MonoInts/README.rst +++ b/src/MonoInts/README.rst @@ -17,96 +17,115 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`ao_mono_elec_integral `_ +`ao_mono_elec_integral `_ array of the mono electronic hamiltonian on the AOs basis : sum of the kinetic and nuclear electronic potential -`ao_overlap `_ +`ao_overlap `_ Overlap between atomic basis functions: :math:`\int \chi_i(r) \chi_j(r) dr)` -`ao_overlap_abs `_ +`ao_overlap_abs `_ Overlap between absolute value of atomic basis functions: :math:`\int |\chi_i(r)| |\chi_j(r)| dr)` -`ao_overlap_x `_ +`ao_overlap_x `_ Overlap between atomic basis functions: :math:`\int \chi_i(r) \chi_j(r) dr)` -`ao_overlap_y `_ +`ao_overlap_y `_ Overlap between atomic basis functions: :math:`\int \chi_i(r) \chi_j(r) dr)` -`ao_overlap_z `_ +`ao_overlap_z `_ Overlap between atomic basis functions: :math:`\int \chi_i(r) \chi_j(r) dr)` -`check_ortho `_ -None -`do_print `_ -None -`n_pt_max_i_x `_ -None -`n_pt_max_integrals `_ -None -`ao_deriv2_x `_ +`check_ortho `_ + Undocumented + +`do_print `_ + Undocumented + +`n_pt_max_i_x `_ + Undocumented + +`n_pt_max_integrals `_ + Undocumented + +`ao_deriv2_x `_ second derivatives matrix elements in the ao basis .. math:: .br {\tt ao_deriv2_x} = \langle \chi_i(x,y,z) \frac{\partial^2}{\partial x^2} |\chi_j (x,y,z) \rangle -`ao_deriv2_y `_ +`ao_deriv2_y `_ second derivatives matrix elements in the ao basis .. math:: .br {\tt ao_deriv2_x} = \langle \chi_i(x,y,z) \frac{\partial^2}{\partial x^2} |\chi_j (x,y,z) \rangle -`ao_deriv2_z `_ +`ao_deriv2_z `_ second derivatives matrix elements in the ao basis .. math:: .br {\tt ao_deriv2_x} = \langle \chi_i(x,y,z) \frac{\partial^2}{\partial x^2} |\chi_j (x,y,z) \rangle -`ao_kinetic_integral `_ +`ao_kinetic_integral `_ array of the priminitve basis kinetic integrals \langle \chi_i |\hat{T}| \chi_j \rangle -`mo_kinetic_integral `_ -None -`mo_mono_elec_integral `_ +`mo_kinetic_integral `_ + Undocumented + +`mo_mono_elec_integral `_ array of the mono electronic hamiltonian on the MOs basis : sum of the kinetic and nuclear electronic potential -`mo_overlap `_ -None -`orthonormalize_mos `_ -None -`ao_nucl_elec_integral `_ +`mo_overlap `_ + Undocumented + +`orthonormalize_mos `_ + Undocumented + +`ao_nucl_elec_integral `_ interaction nuclear electron -`give_polynom_mult_center_mono_elec `_ -None -`i_x1_pol_mult_mono_elec `_ -None -`i_x2_pol_mult_mono_elec `_ -None -`int_gaus_pol `_ -None -`nai_pol_mult `_ -None -`v_e_n `_ -None -`v_phi `_ -None -`v_r `_ -None -`v_theta `_ -None -`wallis `_ -None -`mo_nucl_elec_integral `_ -None -`save_ortho_mos `_ -None +`give_polynom_mult_center_mono_elec `_ + Undocumented + +`i_x1_pol_mult_mono_elec `_ + Undocumented + +`i_x2_pol_mult_mono_elec `_ + Undocumented + +`int_gaus_pol `_ + Undocumented + +`nai_pol_mult `_ + Undocumented + +`v_e_n `_ + Undocumented + +`v_phi `_ + Undocumented + +`v_r `_ + Undocumented + +`v_theta `_ + Undocumented + +`wallis `_ + Undocumented + +`mo_nucl_elec_integral `_ + Undocumented + +`save_ortho_mos `_ + Undocumented + diff --git a/src/MonoInts/spread_dipole_ao.irp.f b/src/MonoInts/spread_dipole_ao.irp.f new file mode 100644 index 00000000..2628c6a0 --- /dev/null +++ b/src/MonoInts/spread_dipole_ao.irp.f @@ -0,0 +1,416 @@ + BEGIN_PROVIDER [ double precision, ao_spread_x, (ao_num_align,ao_num)] + &BEGIN_PROVIDER [ double precision, ao_spread_y, (ao_num_align,ao_num)] + &BEGIN_PROVIDER [ double precision, ao_spread_z, (ao_num_align,ao_num)] + BEGIN_DOC + ! array of the integrals of AO_i * x^2 AO_j + ! array of the integrals of AO_i * y^2 AO_j + ! array of the integrals of AO_i * z^2 AO_j + END_DOC + implicit none + integer :: i,j,n,l + double precision :: f, tmp + integer :: dim1 + double precision :: overlap, overlap_x, overlap_y, overlap_z + double precision :: alpha, beta + double precision :: A_center(3), B_center(3) + integer :: power_A(3), power_B(3) + double precision :: lower_exp_val, dx, c,accu_x,accu_y,accu_z + dim1=500 + lower_exp_val = 40.d0 + ao_spread_x= 0.d0 + ao_spread_y= 0.d0 + ao_spread_z= 0.d0 + !$OMP PARALLEL DO SCHEDULE(GUIDED) & + !$OMP DEFAULT(NONE) & + !$OMP PRIVATE(A_center,B_center,power_A,power_B,& + !$OMP overlap_x,overlap_y, overlap_z, overlap, & + !$OMP alpha, beta,i,j,dx,tmp,c,accu_x,accu_y,accu_z) & + !$OMP SHARED(nucl_coord,ao_power,ao_prim_num, & + !$OMP ao_spread_x,ao_spread_y,ao_spread_z,ao_num,ao_coef_transp,ao_nucl, & + !$OMP ao_expo_transp,dim1,lower_exp_val) + do j=1,ao_num + A_center(1) = nucl_coord( ao_nucl(j), 1 ) + A_center(2) = nucl_coord( ao_nucl(j), 2 ) + A_center(3) = nucl_coord( ao_nucl(j), 3 ) + power_A(1) = ao_power( j, 1 ) + power_A(2) = ao_power( j, 2 ) + power_A(3) = ao_power( j, 3 ) + !DEC$ VECTOR ALIGNED + !DEC$ VECTOR ALWAYS + do i= 1,ao_num + B_center(1) = nucl_coord( ao_nucl(i), 1 ) + B_center(2) = nucl_coord( ao_nucl(i), 2 ) + B_center(3) = nucl_coord( ao_nucl(i), 3 ) + power_B(1) = ao_power( i, 1 ) + power_B(2) = ao_power( i, 2 ) + power_B(3) = ao_power( i, 3 ) + accu_x = 0.d0 + accu_y = 0.d0 + accu_z = 0.d0 + do n = 1,ao_prim_num(j) + alpha = ao_expo_transp(n,j) + !DEC$ VECTOR ALIGNED + do l = 1, ao_prim_num(i) + c = ao_coef_transp(n,j)*ao_coef_transp(l,i) + beta = ao_expo_transp(l,i) + call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1) + call overlap_bourrin_spread(A_center(1),B_center(1),alpha,beta,power_A(1),power_B(1),tmp,lower_exp_val,dx,dim1) + accu_x += c*(tmp*overlap_y*overlap_z) + call overlap_bourrin_spread(A_center(2),B_center(2),alpha,beta,power_A(2),power_B(2),tmp,lower_exp_val,dx,dim1) + accu_y += c*(tmp*overlap_x*overlap_z) + call overlap_bourrin_spread(A_center(3),B_center(3),alpha,beta,power_A(3),power_B(3),tmp,lower_exp_val,dx,dim1) + accu_z += c*(tmp*overlap_y*overlap_x) + enddo + enddo + ao_spread_x(i,j) = accu_x + ao_spread_y(i,j) = accu_y + ao_spread_z(i,j) = accu_z + enddo + enddo + !$OMP END PARALLEL DO + END_PROVIDER + + + + BEGIN_PROVIDER [ double precision, ao_dipole_x, (ao_num_align,ao_num)] + &BEGIN_PROVIDER [ double precision, ao_dipole_y, (ao_num_align,ao_num)] + &BEGIN_PROVIDER [ double precision, ao_dipole_z, (ao_num_align,ao_num)] + BEGIN_DOC + ! array of the integrals of AO_i * x AO_j + ! array of the integrals of AO_i * y AO_j + ! array of the integrals of AO_i * z AO_j + END_DOC + implicit none + integer :: i,j,n,l + double precision :: f, tmp + integer :: dim1 + double precision :: overlap, overlap_x, overlap_y, overlap_z,accu_x,accu_y,accu_z + double precision :: alpha, beta + double precision :: A_center(3), B_center(3) + integer :: power_A(3), power_B(3) + double precision :: lower_exp_val, dx, c + dim1=500 + lower_exp_val = 40.d0 + ao_dipole_x= 0.d0 + ao_dipole_y= 0.d0 + ao_dipole_z= 0.d0 + !$OMP PARALLEL DO SCHEDULE(GUIDED) & + !$OMP DEFAULT(NONE) & + !$OMP PRIVATE(A_center,B_center,power_A,power_B,& + !$OMP overlap_x,overlap_y, overlap_z, overlap, & + !$OMP alpha, beta,i,j,dx,tmp,c,accu_x,accu_y,accu_z) & + !$OMP SHARED(nucl_coord,ao_power,ao_prim_num, & + !$OMP ao_dipole_x,ao_dipole_y,ao_dipole_z,ao_num,ao_coef_transp,ao_nucl, & + !$OMP ao_expo_transp,dim1,lower_exp_val) + do j=1,ao_num + A_center(1) = nucl_coord( ao_nucl(j), 1 ) + A_center(2) = nucl_coord( ao_nucl(j), 2 ) + A_center(3) = nucl_coord( ao_nucl(j), 3 ) + power_A(1) = ao_power( j, 1 ) + power_A(2) = ao_power( j, 2 ) + power_A(3) = ao_power( j, 3 ) + !DEC$ VECTOR ALIGNED + !DEC$ VECTOR ALWAYS + do i= 1,ao_num + B_center(1) = nucl_coord( ao_nucl(i), 1 ) + B_center(2) = nucl_coord( ao_nucl(i), 2 ) + B_center(3) = nucl_coord( ao_nucl(i), 3 ) + power_B(1) = ao_power( i, 1 ) + power_B(2) = ao_power( i, 2 ) + power_B(3) = ao_power( i, 3 ) + accu_x = 0.d0 + accu_y = 0.d0 + accu_z = 0.d0 + do n = 1,ao_prim_num(j) + alpha = ao_expo_transp(n,j) + !DEC$ VECTOR ALIGNED + do l = 1, ao_prim_num(i) + beta = ao_expo_transp(l,i) + c = ao_coef_transp(l,i)*ao_coef_transp(n,j) + call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1) + + call overlap_bourrin_dipole(A_center(1),B_center(1),alpha,beta,power_A(1),power_B(1),tmp,lower_exp_val,dx,dim1) + accu_x = accu_x + c*(tmp*overlap_y*overlap_z) + call overlap_bourrin_dipole(A_center(2),B_center(2),alpha,beta,power_A(2),power_B(2),tmp,lower_exp_val,dx,dim1) + accu_y = accu_y + c*(tmp*overlap_x*overlap_z) + call overlap_bourrin_dipole(A_center(3),B_center(3),alpha,beta,power_A(3),power_B(3),tmp,lower_exp_val,dx,dim1) + accu_z = accu_z + c*(tmp*overlap_y*overlap_x) + enddo + enddo + ao_dipole_x(i,j) = accu_x + ao_dipole_y(i,j) = accu_y + ao_dipole_z(i,j) = accu_z + enddo + enddo + !$OMP END PARALLEL DO + END_PROVIDER + + BEGIN_PROVIDER [ double precision, ao_deriv_1_x, (ao_num_align,ao_num)] + &BEGIN_PROVIDER [ double precision, ao_deriv_1_y, (ao_num_align,ao_num)] + &BEGIN_PROVIDER [ double precision, ao_deriv_1_z, (ao_num_align,ao_num)] + BEGIN_DOC + ! array of the integrals of AO_i * d/dx AO_j + ! array of the integrals of AO_i * d/dy AO_j + ! array of the integrals of AO_i * d/dz AO_j + END_DOC + implicit none + integer :: i,j,n,l + double precision :: f, tmp + integer :: dim1 + double precision :: overlap, overlap_x, overlap_y, overlap_z + double precision :: alpha, beta + double precision :: A_center(3), B_center(3) + integer :: power_A(3), power_B(3) + double precision :: lower_exp_val, dx, c,accu_x,accu_y,accu_z + integer :: i_component + dim1=500 + lower_exp_val = 40.d0 + ao_deriv_1_x= 0.d0 + ao_deriv_1_y= 0.d0 + ao_deriv_1_z= 0.d0 + !$OMP PARALLEL DO SCHEDULE(GUIDED) & + !$OMP DEFAULT(NONE) & + !$OMP PRIVATE(A_center,B_center,power_A,power_B,& + !$OMP overlap_x,overlap_y, overlap_z, overlap, & + !$OMP alpha, beta,i,j,dx,tmp,c,i_component,accu_x,accu_y,accu_z) & + !$OMP SHARED(nucl_coord,ao_power,ao_prim_num, & + !$OMP ao_deriv_1_x,ao_deriv_1_y,ao_deriv_1_z,ao_num,ao_coef_transp,ao_nucl, & + !$OMP ao_expo_transp,dim1,lower_exp_val) + do j=1,ao_num + A_center(1) = nucl_coord( ao_nucl(j), 1 ) + A_center(2) = nucl_coord( ao_nucl(j), 2 ) + A_center(3) = nucl_coord( ao_nucl(j), 3 ) + power_A(1) = ao_power( j, 1 ) + power_A(2) = ao_power( j, 2 ) + power_A(3) = ao_power( j, 3 ) + !DEC$ VECTOR ALIGNED + !DEC$ VECTOR ALWAYS + do i= 1,ao_num + B_center(1) = nucl_coord( ao_nucl(i), 1 ) + B_center(2) = nucl_coord( ao_nucl(i), 2 ) + B_center(3) = nucl_coord( ao_nucl(i), 3 ) + power_B(1) = ao_power( i, 1 ) + power_B(2) = ao_power( i, 2 ) + power_B(3) = ao_power( i, 3 ) + accu_x = 0.d0 + accu_y = 0.d0 + accu_z = 0.d0 + do n = 1,ao_prim_num(j) + alpha = ao_expo_transp(n,j) + !DEC$ VECTOR ALIGNED + do l = 1, ao_prim_num(i) + beta = ao_expo_transp(l,i) + call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1) + c = ao_coef_transp(l,i) * ao_coef_transp(n,j) + i_component = 1 + call overlap_bourrin_deriv_x(i_component,A_center,B_center,alpha,beta,power_A,power_B,dx,lower_exp_val,tmp,dim1) + accu_x += c*(tmp*overlap_y*overlap_z) + i_component = 2 + call overlap_bourrin_deriv_x(i_component,A_center,B_center,alpha,beta,power_A,power_B,dx,lower_exp_val,tmp,dim1) + accu_y += c*(tmp*overlap_x*overlap_z) + i_component = 3 + call overlap_bourrin_deriv_x(i_component,A_center,B_center,alpha,beta,power_A,power_B,dx,lower_exp_val,tmp,dim1) + accu_z += c*(tmp*overlap_y*overlap_x) + enddo + enddo + ao_deriv_1_x(i,j) = accu_x + ao_deriv_1_y(i,j) = accu_y + ao_deriv_1_z(i,j) = accu_z + enddo + enddo + !$OMP END PARALLEL DO + END_PROVIDER + + + + subroutine overlap_bourrin_x_abs(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,lower_exp_val,dx,nx) + implicit none +! compute the following integral : +! int [-infty ; +infty] of [(x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) ] + integer :: i,j,k,l + integer,intent(in) :: power_A,power_B + double precision, intent(in) :: lower_exp_val + double precision,intent(in) :: A_center, B_center,alpha,beta + double precision, intent(out) :: overlap_x,dx + integer, intent(in) :: nx + double precision :: x_min,x_max,domain,x,power,factor,dist,p,p_inv,rho + double precision :: P_center,pouet_timy + if(power_A.lt.0.or.power_B.lt.0)then + overlap_x = 0.d0 + dx = 0.d0 + return + endif + p = alpha + beta + p_inv= 1.d0/p + rho = alpha * beta * p_inv + dist = (A_center - B_center)*(A_center - B_center) + P_center = (alpha * A_center + beta * B_center) * p_inv + factor = dexp(-rho * dist) + pouet_timy = dsqrt(lower_exp_val/p) + x_min = P_center - pouet_timy + x_max = P_center + pouet_timy +!print*,'xmin = ',x_min +!print*,'xmax = ',x_max + domain = x_max-x_min + dx = domain/dble(nx) + overlap_x = 0.d0 + x = x_min + do i = 1, nx + x += dx + overlap_x += (power(power_A,x-A_center) * power(power_B,x-B_center)) * dexp(-p * (x-P_center)*(x-P_center)) + enddo + overlap_x *= factor * dx +end + + subroutine overlap_bourrin_spread(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,lower_exp_val,dx,nx) +! compute the following integral : +! int [-infty ; +infty] of [(x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) * x ] +! needed for the dipole and those things + implicit none + integer :: i,j,k,l + integer,intent(in) :: power_A,power_B + double precision, intent(in) :: lower_exp_val + double precision,intent(in) :: A_center, B_center,alpha,beta + double precision, intent(out) :: overlap_x,dx + integer, intent(in) :: nx + double precision :: x_min,x_max,domain,x,power,factor,dist,p,p_inv,rho + double precision :: P_center,pouet_timy + if(power_A.lt.0.or.power_B.lt.0)then + overlap_x = 0.d0 + dx = 0.d0 + return + endif + p = alpha + beta + p_inv= 1.d0/p + rho = alpha * beta * p_inv + dist = (A_center - B_center)*(A_center - B_center) + P_center = (alpha * A_center + beta * B_center) * p_inv + factor = dexp(-rho * dist) + if(factor.lt.0.000001d0)then +! print*,'factor = ',factor + dx = 0.d0 + overlap_x = 0.d0 + return + endif + pouet_timy = dsqrt(lower_exp_val/p) + x_min = P_center - pouet_timy + x_max = P_center + pouet_timy + domain = x_max-x_min + dx = domain/dble(nx) + overlap_x = 0.d0 + x = x_min + do i = 1, nx + x += dx + overlap_x += power(power_A,x-A_center) * power(power_B,x-B_center) * dexp(-p * (x-P_center)*(x-P_center)) * x * x + enddo + overlap_x *= factor * dx + + end + + double precision function power(n,x) + implicit none + integer :: i,n + double precision :: x,accu + power = x**n + return + end + + subroutine overlap_bourrin_dipole(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,lower_exp_val,dx,nx) +! compute the following integral : +! int [-infty ; +infty] of [(x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) * x ] +! needed for the dipole and those things + implicit none + integer :: i,j,k,l + integer,intent(in) :: power_A,power_B + double precision, intent(in) :: lower_exp_val + double precision,intent(in) :: A_center, B_center,alpha,beta + double precision, intent(out) :: overlap_x,dx + integer, intent(in) :: nx + double precision :: x_min,x_max,domain,x,power,factor,dist,p,p_inv,rho + double precision :: P_center + if(power_A.lt.0.or.power_B.lt.0)then + overlap_x = 0.d0 + dx = 0.d0 + return + endif + p = alpha + beta + p_inv= 1.d0/p + rho = alpha * beta * p_inv + dist = (A_center - B_center)*(A_center - B_center) + P_center = (alpha * A_center + beta * B_center) * p_inv + factor = dexp(-rho * dist) + double precision :: pouet_timy + + pouet_timy = dsqrt(lower_exp_val/p) + x_min = P_center - pouet_timy + x_max = P_center + pouet_timy + domain = x_max-x_min + dx = domain/dble(nx) + overlap_x = 0.d0 + x = x_min + do i = 1, nx + x += dx + overlap_x += power(power_A,x-A_center) * power(power_B,x-B_center) * dexp(-p * (x-P_center)*(x-P_center)) * x + enddo + overlap_x *= factor * dx + + end + + subroutine overlap_bourrin_deriv_x(i_component,A_center,B_center,alpha,beta,power_A,power_B,dx,lower_exp_val,overlap_x,nx) + implicit none + integer :: i,j,k,l + integer,intent(in) :: power_A(3),power_B(3),i_component + double precision,intent(in) :: A_center(3), B_center(3),alpha,beta,lower_exp_val + double precision, intent(out) :: overlap_x,dx + integer, intent(in) :: nx + double precision :: overlap_first, overlap_second +! computes : = (a_x_i - 2 alpha ) + + call overlap_bourrin_x(A_center(i_component),B_center(i_component),alpha,beta,power_A(i_component)-1,power_B(i_component),overlap_first,lower_exp_val,dx,nx) + call overlap_bourrin_x(A_center(i_component),B_center(i_component),alpha,beta,power_A(i_component)+1,power_B(i_component),overlap_second,lower_exp_val,dx,nx) + overlap_x = (power_A(i_component) * overlap_first - 2.d0 * alpha * overlap_second) + end + + subroutine overlap_bourrin_x(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,lower_exp_val,dx,nx) + implicit none +! compute the following integral : +! int [-infty ; +infty] of [(x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) ] + integer :: i,j,k,l + integer,intent(in) :: power_A,power_B + double precision, intent(in) :: lower_exp_val + double precision,intent(in) :: A_center, B_center,alpha,beta + double precision, intent(out) :: overlap_x,dx + integer, intent(in) :: nx + double precision :: x_min,x_max,domain,x,power,factor,dist,p,p_inv,rho + double precision :: P_center,pouet_timy + if(power_A.lt.0.or.power_B.lt.0)then + overlap_x = 0.d0 + dx = 0.d0 + return + endif + p = alpha + beta + p_inv= 1.d0/p + rho = alpha * beta * p_inv + dist = (A_center - B_center)*(A_center - B_center) + P_center = (alpha * A_center + beta * B_center) * p_inv + factor = dexp(-rho * dist) + if(factor.lt.0.000001d0)then + dx = 0.d0 + overlap_x = 0.d0 + return + endif + + pouet_timy = dsqrt(lower_exp_val/p) + x_min = P_center - pouet_timy + x_max = P_center + pouet_timy + domain = x_max-x_min + dx = domain/dble(nx) + overlap_x = 0.d0 + x = x_min + do i = 1, nx + x += dx + overlap_x += power(power_A,x-A_center) * power(power_B,x-B_center) * dexp(-p * (x-P_center)*(x-P_center)) + enddo + overlap_x *= factor * dx + end + diff --git a/src/MonoInts/spread_dipole_mo.irp.f b/src/MonoInts/spread_dipole_mo.irp.f new file mode 100644 index 00000000..d7306727 --- /dev/null +++ b/src/MonoInts/spread_dipole_mo.irp.f @@ -0,0 +1,101 @@ + BEGIN_PROVIDER [double precision, mo_dipole_x , (mo_tot_num_align,mo_tot_num)] +&BEGIN_PROVIDER [double precision, mo_dipole_y , (mo_tot_num_align,mo_tot_num)] +&BEGIN_PROVIDER [double precision, mo_dipole_z , (mo_tot_num_align,mo_tot_num)] + BEGIN_DOC + ! array of the integrals of MO_i * x MO_j + ! array of the integrals of MO_i * y MO_j + ! array of the integrals of MO_i * z MO_j + END_DOC + implicit none + integer :: i1,j1,i,j + double precision :: c_i1,c_j1 + + mo_dipole_x = 0.d0 + mo_dipole_y = 0.d0 + mo_dipole_z = 0.d0 + !$OMP PARALLEL DO DEFAULT(none) & + !$OMP PRIVATE(i,j,i1,j1,c_j1,c_i1) & + !$OMP SHARED(mo_tot_num,ao_num,mo_coef, & + !$OMP mo_dipole_x,mo_dipole_y,mo_dipole_z,ao_dipole_x,ao_dipole_y,ao_dipole_z) + do i = 1, mo_tot_num + do j = 1, mo_tot_num + do i1 = 1,ao_num + c_i1 = mo_coef(i1,i) + do j1 = 1,ao_num + c_j1 = c_i1*mo_coef(j1,j) + mo_dipole_x(j,i) = mo_dipole_x(j,i) + c_j1 * ao_dipole_x(j1,i1) + mo_dipole_y(j,i) = mo_dipole_y(j,i) + c_j1 * ao_dipole_y(j1,i1) + mo_dipole_z(j,i) = mo_dipole_z(j,i) + c_j1 * ao_dipole_z(j1,i1) + enddo + enddo + enddo + enddo + !$OMP END PARALLEL DO +END_PROVIDER + + BEGIN_PROVIDER [double precision, mo_spread_x , (mo_tot_num_align,mo_tot_num)] +&BEGIN_PROVIDER [double precision, mo_spread_y , (mo_tot_num_align,mo_tot_num)] +&BEGIN_PROVIDER [double precision, mo_spread_z , (mo_tot_num_align,mo_tot_num)] + BEGIN_DOC + ! array of the integrals of MO_i * x^2 MO_j + ! array of the integrals of MO_i * y^2 MO_j + ! array of the integrals of MO_i * z^2 MO_j + END_DOC + implicit none + integer :: i1,j1,i,j + double precision :: c_i1,c_j1 + + mo_nucl_elec_integral = 0.d0 + !$OMP PARALLEL DO DEFAULT(none) & + !$OMP PRIVATE(i,j,i1,j1,c_j1,c_i1) & + !$OMP SHARED(mo_tot_num,ao_num,mo_coef, & + !$OMP mo_spread_x,mo_spread_y,mo_spread_z,ao_spread_x,ao_spread_y,ao_spread_z) + do i = 1, mo_tot_num + do j = 1, mo_tot_num + do i1 = 1,ao_num + c_i1 = mo_coef(i1,i) + do j1 = 1,ao_num + c_j1 = c_i1*mo_coef(j1,j) + mo_spread_x(j,i) = mo_spread_x(j,i) + c_j1 * ao_spread_x(j1,i1) + mo_spread_y(j,i) = mo_spread_y(j,i) + c_j1 * ao_spread_y(j1,i1) + mo_spread_z(j,i) = mo_spread_z(j,i) + c_j1 * ao_spread_z(j1,i1) + enddo + enddo + enddo + enddo + !$OMP END PARALLEL DO +END_PROVIDER + + BEGIN_PROVIDER [double precision, mo_deriv_1_x , (mo_tot_num_align,mo_tot_num)] +&BEGIN_PROVIDER [double precision, mo_deriv_1_y , (mo_tot_num_align,mo_tot_num)] +&BEGIN_PROVIDER [double precision, mo_deriv_1_z , (mo_tot_num_align,mo_tot_num)] + BEGIN_DOC + ! array of the integrals of MO_i * d/dx MO_j + ! array of the integrals of MO_i * d/dy MO_j + ! array of the integrals of MO_i * d/dz MO_j + END_DOC + implicit none + integer :: i1,j1,i,j + double precision :: c_i1,c_j1 + + mo_nucl_elec_integral = 0.d0 + !$OMP PARALLEL DO DEFAULT(none) & + !$OMP PRIVATE(i,j,i1,j1,c_j1,c_i1) & + !$OMP SHARED(mo_tot_num,ao_num,mo_coef, & + !$OMP mo_deriv_1_x,mo_deriv_1_y,mo_deriv_1_z,ao_spread_x,ao_spread_y,ao_spread_z) + do i = 1, mo_tot_num + do j = 1, mo_tot_num + do i1 = 1,ao_num + c_i1 = mo_coef(i1,i) + do j1 = 1,ao_num + c_j1 = c_i1*mo_coef(j1,j) + mo_deriv_1_x(j,i) = mo_deriv_1_x(j,i) + c_j1 * ao_spread_x(j1,i1) + mo_deriv_1_y(j,i) = mo_deriv_1_y(j,i) + c_j1 * ao_spread_y(j1,i1) + mo_deriv_1_z(j,i) = mo_deriv_1_z(j,i) + c_j1 * ao_spread_z(j1,i1) + enddo + enddo + enddo + enddo + !$OMP END PARALLEL DO +END_PROVIDER + diff --git a/src/README.rst b/src/README.rst index 9c6b5bcb..3aae2807 100644 --- a/src/README.rst +++ b/src/README.rst @@ -129,4 +129,13 @@ Needed Modules * `Dets `_ * `DensityMatrix `_ * `CISD `_ +* `Perturbation `_ + +Documentation +============= + +.. Do not edit this section. It was auto-generated from the +.. NEEDED_MODULES file. + + diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f index ecae8cdc..0e47233c 100644 --- a/src/Utils/LinearAlgebra.irp.f +++ b/src/Utils/LinearAlgebra.irp.f @@ -181,6 +181,13 @@ subroutine lapack_diag(eigvalues,eigvectors,H,nmax,n) implicit none BEGIN_DOC ! Diagonalize matrix H + ! + ! H is untouched between input and ouptut + ! + ! eigevalues(i) = ith lowest eigenvalue of the H matrix + ! + ! eigvectors(i,j) = where i is the basis function and psi_j is the j th eigenvector + ! END_DOC integer, intent(in) :: n,nmax double precision, intent(out) :: eigvectors(nmax,n) @@ -189,7 +196,6 @@ subroutine lapack_diag(eigvalues,eigvectors,H,nmax,n) double precision,allocatable :: eigenvalues(:) double precision,allocatable :: work(:) double precision,allocatable :: A(:,:) - !eigvectors(i,j) = where d_i is the basis function and psi_j is the j th eigenvector allocate(A(nmax,n),eigenvalues(nmax),work(4*nmax)) integer :: LWORK, info, i,j,l,k A=H