From 6e91ca9104753fde1b1faf007926480a3d0961d6 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Mon, 27 Mar 2017 15:28:20 +0200 Subject: [PATCH] minor modifs mend --- plugins/MRPT/MRPT_Utils.main.irp.f | 15 +-------- plugins/MRPT_Utils/MRMP2_density.irp.f | 46 ++++++++++++++++++++++++++ plugins/MRPT_Utils/mrpt_dress.irp.f | 4 ++- plugins/MRPT_Utils/mrpt_utils.irp.f | 32 +++++++++--------- plugins/Psiref_CAS/psi_ref.irp.f | 14 ++++++-- src/MO_Basis/rotate_mos.irp.f | 9 +++++ src/Utils/invert.irp.f | 19 +++++++++++ 7 files changed, 107 insertions(+), 32 deletions(-) create mode 100644 plugins/MRPT_Utils/MRMP2_density.irp.f create mode 100644 src/MO_Basis/rotate_mos.irp.f create mode 100644 src/Utils/invert.irp.f diff --git a/plugins/MRPT/MRPT_Utils.main.irp.f b/plugins/MRPT/MRPT_Utils.main.irp.f index 24361312..1b6efb4f 100644 --- a/plugins/MRPT/MRPT_Utils.main.irp.f +++ b/plugins/MRPT/MRPT_Utils.main.irp.f @@ -46,19 +46,6 @@ end subroutine routine_2 implicit none - integer :: i - do i = 1, n_core_inact_orb - print*,fock_core_inactive_total(i,1,1),fock_core_inactive(i) - enddo - double precision :: accu - accu = 0.d0 - do i = 1, n_act_orb - integer :: j_act_orb - j_act_orb = list_act(i) - accu += one_body_dm_mo_alpha(j_act_orb,j_act_orb,1) - print*,one_body_dm_mo_alpha(j_act_orb,j_act_orb,1),one_body_dm_mo_beta(j_act_orb,j_act_orb,1) - enddo - print*,'accu = ',accu - + provide electronic_psi_ref_average_value end diff --git a/plugins/MRPT_Utils/MRMP2_density.irp.f b/plugins/MRPT_Utils/MRMP2_density.irp.f new file mode 100644 index 00000000..1051edf9 --- /dev/null +++ b/plugins/MRPT_Utils/MRMP2_density.irp.f @@ -0,0 +1,46 @@ +BEGIN_PROVIDER [double precision, MRMP2_density, (mo_tot_num_align, mo_tot_num)] + implicit none + integer :: i,j,k,l + double precision :: accu, mp2_dm(mo_tot_num) + MRMP2_density = one_body_dm_mo + call give_2h2p_density(mp2_dm) + accu = 0.d0 + do i = 1, n_virt_orb + j = list_virt(i) + accu += mp2_dm(j) + MRMP2_density(j,j)+= mp2_dm(j) + enddo + +END_PROVIDER + +subroutine give_2h2p_density(mp2_density_diag_alpha_beta) + implicit none + double precision, intent(out) :: mp2_density_diag_alpha_beta(mo_tot_num) + integer :: i,j,k,l,m + integer :: iorb,jorb,korb,lorb + + double precision :: get_mo_bielec_integral + double precision :: direct_int + double precision :: coef_double + + mp2_density_diag_alpha_beta = 0.d0 + do k = 1, n_virt_orb + korb = list_virt(k) + do i = 1, n_inact_orb + iorb = list_inact(i) + do j = 1, n_inact_orb + jorb = list_inact(j) + do l = 1, n_virt_orb + lorb = list_virt(l) + direct_int = get_mo_bielec_integral(iorb,jorb,korb,lorb ,mo_integrals_map) + coef_double = direct_int/(fock_core_inactive_total_spin_trace(iorb,1) + fock_core_inactive_total_spin_trace(jorb,1) & + -fock_virt_total_spin_trace(korb,1) - fock_virt_total_spin_trace(lorb,1)) + mp2_density_diag_alpha_beta(korb) += coef_double * coef_double + enddo + enddo + enddo + print*, mp2_density_diag_alpha_beta(korb) + enddo + +end + diff --git a/plugins/MRPT_Utils/mrpt_dress.irp.f b/plugins/MRPT_Utils/mrpt_dress.irp.f index 9699a1df..f722ce4c 100644 --- a/plugins/MRPT_Utils/mrpt_dress.irp.f +++ b/plugins/MRPT_Utils/mrpt_dress.irp.f @@ -121,7 +121,8 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip delta_e(i_state) = 1.d+20 enddo else - call get_delta_e_dyall(psi_ref(1,1,index_i),tq(1,1,i_alpha),coef_array,hialpha,delta_e) + call get_delta_e_dyall(psi_ref(1,1,index_i),tq(1,1,i_alpha),delta_e) +! print*, 'delta_e',delta_e !!!!!!!!!!!!! SHIFTED BK ! double precision :: hjj ! call i_h_j(tq(1,1,i_alpha),tq(1,1,i_alpha),Nint,hjj) @@ -129,6 +130,7 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! endif hij_array(index_i) = hialpha +! print*, 'hialpha ',hialpha do i_state = 1,N_states delta_e_inv_array(index_i,i_state) = 1.d0/delta_e(i_state) enddo diff --git a/plugins/MRPT_Utils/mrpt_utils.irp.f b/plugins/MRPT_Utils/mrpt_utils.irp.f index 5ef9516b..c5418847 100644 --- a/plugins/MRPT_Utils/mrpt_utils.irp.f +++ b/plugins/MRPT_Utils/mrpt_utils.irp.f @@ -40,7 +40,6 @@ enddo print*, '1h = ',accu - stop ! 1p delta_ij_tmp = 0.d0 call H_apply_mrpt_1p(delta_ij_tmp,N_det) @@ -235,7 +234,7 @@ END_PROVIDER enddo END_PROVIDER - BEGIN_PROVIDER [ double precision, CI_electronic_dressed_pt2_new_energy, (N_states_diag) ] + BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_electronic_energy, (N_states_diag) ] &BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_eigenvectors, (N_det,N_states_diag) ] &BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_eigenvectors_s2, (N_states_diag) ] BEGIN_DOC @@ -279,7 +278,7 @@ END_PROVIDER allocate (eigenvectors(size(H_matrix_all_dets,1),N_det)) allocate (eigenvalues(N_det)) call lapack_diag(eigenvalues,eigenvectors, & - H_matrix_all_dets,size(H_matrix_all_dets,1),N_det) + Hmatrix_dressed_pt2_new_symmetrized,size(H_matrix_all_dets,1),N_det) CI_electronic_energy(:) = 0.d0 if (s2_eig) then i_state = 0 @@ -288,8 +287,11 @@ END_PROVIDER good_state_array = .False. call u_0_S2_u_0(s2_eigvalues,eigenvectors,N_det,psi_det,N_int,& N_det,size(eigenvectors,1)) +! print*, s2_eigvalues(:) + print*,'N_det',N_det do j=1,N_det ! Select at least n_states states with S^2 values closed to "expected_s2" + print*, s2_eigvalues(j),expected_s2 if(dabs(s2_eigvalues(j)-expected_s2).le.0.5d0)then i_state +=1 index_good_state_array(i_state) = j @@ -303,10 +305,10 @@ END_PROVIDER ! Fill the first "i_state" states that have a correct S^2 value do j = 1, i_state do i=1,N_det - CI_eigenvectors(i,j) = eigenvectors(i,index_good_state_array(j)) + CI_dressed_pt2_new_eigenvectors(i,j) = eigenvectors(i,index_good_state_array(j)) enddo - CI_electronic_energy(j) = eigenvalues(index_good_state_array(j)) - CI_eigenvectors_s2(j) = s2_eigvalues(index_good_state_array(j)) + CI_dressed_pt2_new_electronic_energy(j) = eigenvalues(index_good_state_array(j)) + CI_dressed_pt2_new_eigenvectors_s2(j) = s2_eigvalues(index_good_state_array(j)) enddo i_other_state = 0 do j = 1, N_det @@ -316,10 +318,10 @@ END_PROVIDER exit endif do i=1,N_det - CI_eigenvectors(i,i_state+i_other_state) = eigenvectors(i,j) + CI_dressed_pt2_new_eigenvectors(i,i_state+i_other_state) = eigenvectors(i,j) enddo - CI_electronic_energy(i_state+i_other_state) = eigenvalues(j) - CI_eigenvectors_s2(i_state+i_other_state) = s2_eigvalues(i_state+i_other_state) + CI_dressed_pt2_new_electronic_energy(i_state+i_other_state) = eigenvalues(j) + CI_dressed_pt2_new_eigenvectors_s2(i_state+i_other_state) = s2_eigvalues(i_state+i_other_state) enddo else @@ -334,10 +336,10 @@ END_PROVIDER print*,'' do j=1,min(N_states_diag,N_det) do i=1,N_det - CI_eigenvectors(i,j) = eigenvectors(i,j) + CI_dressed_pt2_new_eigenvectors(i,j) = eigenvectors(i,j) enddo - CI_electronic_energy(j) = eigenvalues(j) - CI_eigenvectors_s2(j) = s2_eigvalues(j) + CI_dressed_pt2_new_electronic_energy(j) = eigenvalues(j) + CI_dressed_pt2_new_eigenvectors_s2(j) = s2_eigvalues(j) enddo endif deallocate(index_good_state_array,good_state_array) @@ -348,9 +350,9 @@ END_PROVIDER ! Select the "N_states_diag" states of lowest energy do j=1,min(N_det,N_states_diag) do i=1,N_det - CI_eigenvectors(i,j) = eigenvectors(i,j) + CI_dressed_pt2_new_eigenvectors(i,j) = eigenvectors(i,j) enddo - CI_electronic_energy(j) = eigenvalues(j) + CI_dressed_pt2_new_electronic_energy(j) = eigenvalues(j) enddo endif deallocate(eigenvectors,eigenvalues) @@ -370,7 +372,7 @@ BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_energy, (N_states_diag) ] character*(8) :: st call write_time(output_determinants) do j=1,N_states_diag - CI_dressed_pt2_new_energy(j) = CI_electronic_dressed_pt2_new_energy(j) + nuclear_repulsion + CI_dressed_pt2_new_energy(j) = CI_dressed_pt2_new_electronic_energy(j) + nuclear_repulsion write(st,'(I4)') j call write_double(output_determinants,CI_dressed_pt2_new_energy(j),'Energy of state '//trim(st)) call write_double(output_determinants,CI_eigenvectors_s2(j),'S^2 of state '//trim(st)) diff --git a/plugins/Psiref_CAS/psi_ref.irp.f b/plugins/Psiref_CAS/psi_ref.irp.f index 0729a540..1f337532 100644 --- a/plugins/Psiref_CAS/psi_ref.irp.f +++ b/plugins/Psiref_CAS/psi_ref.irp.f @@ -72,9 +72,19 @@ END_PROVIDER &BEGIN_PROVIDER [double precision, psi_ref_average_value, (N_states)] implicit none integer :: i,j - call u_0_H_u_0(electronic_psi_ref_average_value,psi_ref_coef,N_det_ref,psi_ref,N_int,N_states,psi_det_size) + electronic_psi_ref_average_value = psi_energy do i = 1, N_states - psi_ref_average_value(i) = electronic_psi_ref_average_value(i) + nuclear_repulsion + psi_ref_average_value(i) = psi_energy(i) + nuclear_repulsion enddo + double precision :: accu,hij + accu = 0.d0 + do i = 1, N_det_ref + do j = 1, N_det_ref + call i_H_j(psi_ref(1,1,i),psi_ref(1,1,j),N_int,hij) + accu += psi_ref_coef(i,1) * psi_ref_coef(j,1) * hij + enddo + enddo + electronic_psi_ref_average_value(1) = accu + psi_ref_average_value(1) = electronic_psi_ref_average_value(1) + nuclear_repulsion END_PROVIDER diff --git a/src/MO_Basis/rotate_mos.irp.f b/src/MO_Basis/rotate_mos.irp.f new file mode 100644 index 00000000..fe664018 --- /dev/null +++ b/src/MO_Basis/rotate_mos.irp.f @@ -0,0 +1,9 @@ +program rotate_mos + implicit none + integer :: i,j + write(*,*)'Which couple of MOs would you like to mix ?' + read(5,*)i,j + call mix_mo_jk(i,j) + call save_mos + +end diff --git a/src/Utils/invert.irp.f b/src/Utils/invert.irp.f new file mode 100644 index 00000000..4c626cca --- /dev/null +++ b/src/Utils/invert.irp.f @@ -0,0 +1,19 @@ +subroutine invert_matrix(A,LDA,na,A_inv,LDA_inv) +implicit none +double precision, intent(in) :: A (LDA,na) +integer, intent(in) :: LDA, LDA_inv +integer, intent(in) :: na +double precision, intent(out) :: A_inv (LDA_inv,na) + + double precision :: work(LDA_inv*max(na,64)) +!DIR$ ATTRIBUTES ALIGN: $IRP_ALIGN :: work + integer :: inf + integer :: ipiv(LDA_inv) +!DIR$ ATTRIBUTES ALIGN: $IRP_ALIGN :: ipiv + integer :: lwork + A_inv(1:na,1:na) = A(1:na,1:na) + call dgetrf(na, na, A_inv, LDA_inv, ipiv, inf ) + lwork = SIZE(work) + call dgetri(na, A_inv, LDA_inv, ipiv, work, lwork, inf ) +end +