diff --git a/plugins/DFT_Utils/grid_density.irp.f b/plugins/DFT_Utils/grid_density.irp.f index 6071a18b..e1296ce0 100644 --- a/plugins/DFT_Utils/grid_density.irp.f +++ b/plugins/DFT_Utils/grid_density.irp.f @@ -5,7 +5,7 @@ END_PROVIDER BEGIN_PROVIDER [integer, n_points_radial_grid] implicit none - n_points_radial_grid = 10000 + n_points_radial_grid = 10 END_PROVIDER @@ -22,21 +22,21 @@ END_PROVIDER integer :: i double precision :: accu double precision :: degre_rad -!degre_rad = 180.d0/pi -!accu = 0.d0 -!do i = 1, n_points_integration_angular_lebedev -! accu += weights_angular_integration_lebedev(i) -! weights_angular_points(i) = weights_angular_integration_lebedev(i) * 2.d0 * pi -! angular_quadrature_points(i,1) = dcos ( degre_rad * theta_angular_integration_lebedev(i)) & -! * dsin ( degre_rad * phi_angular_integration_lebedev(i)) -! angular_quadrature_points(i,2) = dsin ( degre_rad * theta_angular_integration_lebedev(i)) & -! * dsin ( degre_rad * phi_angular_integration_lebedev(i)) -! angular_quadrature_points(i,3) = dcos ( degre_rad * phi_angular_integration_lebedev(i)) -!enddo -!print*,'ANGULAR' -!print*,'' -!print*,'accu = ',accu -!ASSERT( dabs(accu - 1.D0) < 1.d-10) + degre_rad = 180.d0/pi + accu = 0.d0 + do i = 1, n_points_integration_angular_lebedev + accu += weights_angular_integration_lebedev(i) + weights_angular_points(i) = weights_angular_integration_lebedev(i) * 2.d0 * pi + angular_quadrature_points(i,1) = dcos ( degre_rad * theta_angular_integration_lebedev(i)) & + * dsin ( degre_rad * phi_angular_integration_lebedev(i)) + angular_quadrature_points(i,2) = dsin ( degre_rad * theta_angular_integration_lebedev(i)) & + * dsin ( degre_rad * phi_angular_integration_lebedev(i)) + angular_quadrature_points(i,3) = dcos ( degre_rad * phi_angular_integration_lebedev(i)) + enddo + print*,'ANGULAR' + print*,'' + print*,'accu = ',accu + ASSERT( dabs(accu - 1.D0) < 1.d-10) END_PROVIDER @@ -152,8 +152,8 @@ END_PROVIDER do i = 1, mo_tot_num do m = 1, mo_tot_num contrib = mos_array(i) * mos_array(m) - one_body_dm_mo_alpha_at_grid_points(l,k,j) += one_body_dm_mo_alpha(i,m) * contrib - one_body_dm_mo_beta_at_grid_points(l,k,j) += one_body_dm_mo_beta(i,m) * contrib + one_body_dm_mo_alpha_at_grid_points(l,k,j) += one_body_dm_mo_alpha_average(i,m) * contrib + one_body_dm_mo_beta_at_grid_points(l,k,j) += one_body_dm_mo_beta_average(i,m) * contrib enddo enddo diff --git a/plugins/MRPT_Utils/energies_cas.irp.f b/plugins/MRPT_Utils/energies_cas.irp.f index 02077ebc..ef026e07 100644 --- a/plugins/MRPT_Utils/energies_cas.irp.f +++ b/plugins/MRPT_Utils/energies_cas.irp.f @@ -293,7 +293,8 @@ BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2 END_PROVIDER -BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_act_orb,2,2,2,N_states)] + BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_act_orb,2,2,2,N_states)] +&BEGIN_PROVIDER [ double precision, two_anhil_one_creat_norm, (n_act_orb,n_act_orb,n_act_orb,2,2,2,N_states)] implicit none integer :: i,j integer :: ispin,jspin,kspin @@ -344,6 +345,7 @@ BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_a norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states) call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) two_anhil_one_creat(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall_no_exchange(state_target) - energies(state_target) + two_anhil_one_creat_norm(iorb,jorb,korb,ispin,jspin,kspin,state_target) = norm_out(state_target) enddo enddo enddo @@ -355,7 +357,54 @@ BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_a END_PROVIDER -BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_act_orb,2,2,2,N_states)] + + BEGIN_PROVIDER [ double precision, two_anhil_one_creat_spin_average, (n_act_orb,n_act_orb,n_act_orb,N_states)] + implicit none + integer :: i,j + integer :: ispin,jspin,kspin + integer :: orb_i, hole_particle_i,spin_exc_i + integer :: orb_j, hole_particle_j,spin_exc_j + integer :: orb_k, hole_particle_k,spin_exc_k + double precision :: norm_out(N_states) + integer(bit_kind), allocatable :: psi_in_out(:,:,:) + double precision, allocatable :: psi_in_out_coef(:,:) + use bitmasks + allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det_ref,N_states)) + + integer :: iorb,jorb + integer :: korb + integer :: state_target + double precision :: energies(n_states) + double precision :: accu + do iorb = 1,n_act_orb + orb_i = list_act(iorb) + do jorb = 1, n_act_orb + orb_j = list_act(jorb) + do korb = 1, n_act_orb + orb_k = list_act(korb) + do state_target = 1, N_states + accu = 0.d0 + do ispin = 1,2 + do jspin = 1,2 + do kspin = 1,2 + two_anhil_one_creat_spin_average(iorb,jorb,korb,state_target) += two_anhil_one_creat(iorb,jorb,korb,ispin,jspin,kspin,state_target)* & + two_anhil_one_creat_norm(iorb,jorb,korb,ispin,jspin,kspin,state_target) + accu += two_anhil_one_creat_norm(iorb,jorb,korb,ispin,jspin,kspin,state_target) + enddo + enddo + enddo + two_anhil_one_creat_spin_average(iorb,jorb,korb,state_target) = two_anhil_one_creat_spin_average(iorb,jorb,korb,state_target) /accu + enddo + enddo + enddo + enddo + deallocate(psi_in_out,psi_in_out_coef) + +END_PROVIDER + + + BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_act_orb,2,2,2,N_states)] +&BEGIN_PROVIDER [ double precision, two_creat_one_anhil_norm, (n_act_orb,n_act_orb,n_act_orb,2,2,2,N_states)] implicit none integer :: i,j integer :: ispin,jspin,kspin @@ -406,6 +455,8 @@ implicit none norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states) call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) two_creat_one_anhil(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall_no_exchange(state_target) - energies(state_target) + two_creat_one_anhil_norm(iorb,jorb,korb,ispin,jspin,kspin,state_target) = norm_out(state_target) +! print*, norm_out(state_target) enddo enddo enddo @@ -415,6 +466,51 @@ implicit none enddo deallocate(psi_in_out,psi_in_out_coef) +END_PROVIDER + + + BEGIN_PROVIDER [ double precision, two_creat_one_anhil_spin_average, (n_act_orb,n_act_orb,n_act_orb,N_states)] +implicit none + integer :: i,j + integer :: ispin,jspin,kspin + integer :: orb_i, hole_particle_i,spin_exc_i + integer :: orb_j, hole_particle_j,spin_exc_j + integer :: orb_k, hole_particle_k,spin_exc_k + double precision :: norm_out(N_states) + integer(bit_kind), allocatable :: psi_in_out(:,:,:) + double precision, allocatable :: psi_in_out_coef(:,:) + use bitmasks + allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det_ref,N_states)) + + integer :: iorb,jorb + integer :: korb + integer :: state_target + double precision :: energies(n_states),accu + do iorb = 1,n_act_orb + orb_i = list_act(iorb) + do jorb = 1, n_act_orb + orb_j = list_act(jorb) + do korb = 1, n_act_orb + orb_k = list_act(korb) + do state_target = 1, N_states + accu = 0.d0 + do ispin = 1,2 + do jspin = 1,2 + do kspin = 1,2 + two_creat_one_anhil_spin_average(iorb,jorb,korb,state_target) += two_creat_one_anhil(iorb,jorb,korb,ispin,jspin,kspin,state_target) * & + two_creat_one_anhil_norm(iorb,jorb,korb,ispin,jspin,kspin,state_target) + accu += two_creat_one_anhil_norm(iorb,jorb,korb,ispin,jspin,kspin,state_target) + print*, accu + enddo + enddo + enddo + two_creat_one_anhil_spin_average(iorb,jorb,korb,state_target) = two_creat_one_anhil_spin_average(iorb,jorb,korb,state_target) / accu + enddo + enddo + enddo + enddo + deallocate(psi_in_out,psi_in_out_coef) + END_PROVIDER !BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_act_orb,N_states)] diff --git a/plugins/MRPT_Utils/mrpt_dress.irp.f b/plugins/MRPT_Utils/mrpt_dress.irp.f index f722ce4c..a08b6108 100644 --- a/plugins/MRPT_Utils/mrpt_dress.irp.f +++ b/plugins/MRPT_Utils/mrpt_dress.irp.f @@ -122,11 +122,14 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip enddo else call get_delta_e_dyall(psi_ref(1,1,index_i),tq(1,1,i_alpha),delta_e) + if(degree_scalar.eq.1)then + delta_e = 1.d+20 + endif ! 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) -! delta_e(1) = CI_electronic_energy(1) - hjj +! delta_e(1) = electronic_psi_ref_average_value(1) - hjj ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! endif hij_array(index_i) = hialpha diff --git a/plugins/MRPT_Utils/mrpt_utils.irp.f b/plugins/MRPT_Utils/mrpt_utils.irp.f index c5418847..79aa624f 100644 --- a/plugins/MRPT_Utils/mrpt_utils.irp.f +++ b/plugins/MRPT_Utils/mrpt_utils.irp.f @@ -40,40 +40,38 @@ enddo print*, '1h = ',accu - ! 1p - delta_ij_tmp = 0.d0 - call H_apply_mrpt_1p(delta_ij_tmp,N_det) - accu = 0.d0 - do i_state = 1, N_states - do i = 1, N_det - do j = 1, N_det - accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) - delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) - enddo - write(*,'(1000(F16.10,x))')delta_ij_tmp(i,:,:) - enddo - second_order_pt_new_1p(i_state) = accu(i_state) - enddo - print*, '1p = ',accu +!! 1p +!delta_ij_tmp = 0.d0 +!call H_apply_mrpt_1p(delta_ij_tmp,N_det) +!accu = 0.d0 +!do i_state = 1, N_states +!do i = 1, N_det +! do j = 1, N_det +! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) +! delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) +! enddo +! write(*,'(1000(F16.10,x))')delta_ij_tmp(i,:,:) +!enddo +!second_order_pt_new_1p(i_state) = accu(i_state) +!enddo +!print*, '1p = ',accu ! 1h1p - delta_ij_tmp = 0.d0 - call H_apply_mrpt_1h1p(delta_ij_tmp,N_det) - double precision :: e_corr_from_1h1p_singles(N_states) -!call give_singles_and_partial_doubles_1h1p_contrib(delta_ij_tmp,e_corr_from_1h1p_singles) -!call give_1h1p_only_doubles_spin_cross(delta_ij_tmp) - accu = 0.d0 - do i_state = 1, N_states - do i = 1, N_det - do j = 1, N_det - accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) - delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) - enddo - write(*,'(1000(F16.10,x))')delta_ij_tmp(i,:,:) - enddo - second_order_pt_new_1h1p(i_state) = accu(i_state) - enddo - print*, '1h1p = ',accu +!delta_ij_tmp = 0.d0 +!call H_apply_mrpt_1h1p(delta_ij_tmp,N_det) +!double precision :: e_corr_from_1h1p_singles(N_states) +!accu = 0.d0 +!do i_state = 1, N_states +!do i = 1, N_det +! do j = 1, N_det +! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) +! delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) +! enddo +! write(*,'(1000(F16.10,x))')delta_ij_tmp(i,:,:) +!enddo +!second_order_pt_new_1h1p(i_state) = accu(i_state) +!enddo +!print*, '1h1p = ',accu ! 1h1p third order if(do_third_order_1h1p)then @@ -93,73 +91,73 @@ print*, '1h1p(3)',accu endif - ! 2h - delta_ij_tmp = 0.d0 - call H_apply_mrpt_2h(delta_ij_tmp,N_det) - accu = 0.d0 - do i_state = 1, N_states - do i = 1, N_det - do j = 1, N_det - accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) - delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) - enddo - write(*,'(1000(F16.10,x))')delta_ij_tmp(i,:,:) - enddo - second_order_pt_new_2h(i_state) = accu(i_state) - enddo - print*, '2h = ',accu +!! 2h +!delta_ij_tmp = 0.d0 +!call H_apply_mrpt_2h(delta_ij_tmp,N_det) +!accu = 0.d0 +!do i_state = 1, N_states +!do i = 1, N_det +! do j = 1, N_det +! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) +! delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) +! enddo +! write(*,'(1000(F16.10,x))')delta_ij_tmp(i,:,:) +!enddo +!second_order_pt_new_2h(i_state) = accu(i_state) +!enddo +!print*, '2h = ',accu - ! 2p - delta_ij_tmp = 0.d0 - call H_apply_mrpt_2p(delta_ij_tmp,N_det) - accu = 0.d0 - do i_state = 1, N_states - do i = 1, N_det - do j = 1, N_det - accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) - delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) - enddo - write(*,'(1000(F16.10,x))')delta_ij_tmp(i,:,:) - enddo - second_order_pt_new_2p(i_state) = accu(i_state) - enddo - print*, '2p = ',accu +!! 2p +!delta_ij_tmp = 0.d0 +!call H_apply_mrpt_2p(delta_ij_tmp,N_det) +!accu = 0.d0 +!do i_state = 1, N_states +!do i = 1, N_det +! do j = 1, N_det +! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) +! delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) +! enddo +! write(*,'(1000(F16.10,x))')delta_ij_tmp(i,:,:) +!enddo +!second_order_pt_new_2p(i_state) = accu(i_state) +!enddo +!print*, '2p = ',accu ! 1h2p delta_ij_tmp = 0.d0 !call give_1h2p_contrib(delta_ij_tmp) - call H_apply_mrpt_1h2p(delta_ij_tmp,N_det) - accu = 0.d0 - do i_state = 1, N_states - do i = 1, N_det - do j = 1, N_det - accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) - delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) - enddo - write(*,'(1000(F16.10,x))')delta_ij_tmp(i,:,:) - enddo - second_order_pt_new_1h2p(i_state) = accu(i_state) - enddo - print*, '1h2p = ',accu +!call H_apply_mrpt_1h2p(delta_ij_tmp,N_det) +!accu = 0.d0 +!do i_state = 1, N_states +!do i = 1, N_det +! do j = 1, N_det +! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) +! delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) +! enddo +! write(*,'(1000(F16.10,x))')delta_ij_tmp(i,:,:) +!enddo +!second_order_pt_new_1h2p(i_state) = accu(i_state) +!enddo +!print*, '1h2p = ',accu - ! 2h1p - delta_ij_tmp = 0.d0 +!! 2h1p +!delta_ij_tmp = 0.d0 !call give_2h1p_contrib(delta_ij_tmp) - call H_apply_mrpt_2h1p(delta_ij_tmp,N_det) - accu = 0.d0 - do i_state = 1, N_states - do i = 1, N_det - do j = 1, N_det - accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) - delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) - enddo - write(*,'(1000(F16.10,x))')delta_ij_tmp(i,:,:) - enddo - second_order_pt_new_2h1p(i_state) = accu(i_state) - enddo - print*, '2h1p = ',accu +!call H_apply_mrpt_2h1p(delta_ij_tmp,N_det) +!accu = 0.d0 +!do i_state = 1, N_states +!do i = 1, N_det +! do j = 1, N_det +! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) +! delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) +! enddo +! write(*,'(1000(F16.10,x))')delta_ij_tmp(i,:,:) +!enddo +!second_order_pt_new_2h1p(i_state) = accu(i_state) +!enddo +!print*, '2h1p = ',accu - ! 2h2p +!! 2h2p !delta_ij_tmp = 0.d0 !call H_apply_mrpt_2h2p(delta_ij_tmp,N_det) !accu = 0.d0 @@ -287,7 +285,6 @@ 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" diff --git a/plugins/MRPT_Utils/psi_active_prov.irp.f b/plugins/MRPT_Utils/psi_active_prov.irp.f index d9c5fda3..f86947d8 100644 --- a/plugins/MRPT_Utils/psi_active_prov.irp.f +++ b/plugins/MRPT_Utils/psi_active_prov.irp.f @@ -354,7 +354,8 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final) kspin = particle_list_practical(1,1) i_particle_act = particle_list_practical(2,1) do i_state = 1, N_states - delta_e_act(i_state) += two_anhil_one_creat(i_particle_act,i_hole_act,j_hole_act,kspin,ispin,jspin,i_state) +! delta_e_act(i_state) += two_anhil_one_creat(i_particle_act,i_hole_act,j_hole_act,kspin,ispin,jspin,i_state) + delta_e_act(i_state) += two_anhil_one_creat_spin_average(i_particle_act,i_hole_act,j_hole_act,i_state) enddo else if (n_holes_act == 1 .and. n_particles_act == 2) then @@ -369,7 +370,9 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final) j_particle_act = particle_list_practical(2,2) do i_state = 1, N_states - delta_e_act(i_state) += two_creat_one_anhil(i_particle_act,j_particle_act,i_hole_act,jspin,kspin,ispin,i_state) +! delta_e_act(i_state) += two_creat_one_anhil(i_particle_act,j_particle_act,i_hole_act,jspin,kspin,ispin,i_state) + delta_e_act(i_state) += 0.5d0 * (two_creat_one_anhil_spin_average(i_particle_act,j_particle_act,i_hole_act,i_state) & + +two_creat_one_anhil_spin_average(j_particle_act,i_particle_act,i_hole_act,i_state)) enddo else if (n_holes_act == 3 .and. n_particles_act == 0) then diff --git a/src/MO_Basis/cholesky_mo.irp.f b/src/MO_Basis/cholesky_mo.irp.f index 5d29db41..774198a3 100644 --- a/src/MO_Basis/cholesky_mo.irp.f +++ b/src/MO_Basis/cholesky_mo.irp.f @@ -50,9 +50,9 @@ subroutine cholesky_mo(n,m,P,LDP,C,LDC,tol_in,rank) deallocate(W,work) end -subroutine svd_mo(n,m,P,LDP,C,LDC) - implicit none - BEGIN_DOC +!subroutine svd_mo(n,m,P,LDP,C,LDC) +!implicit none +!BEGIN_DOC ! Singular value decomposition of the AO Density matrix ! ! n : Number of AOs @@ -67,22 +67,22 @@ subroutine svd_mo(n,m,P,LDP,C,LDC) ! ! rank : Nomber of local MOs (output) ! - END_DOC - integer, intent(in) :: n,m, LDC, LDP - double precision, intent(in) :: P(LDP,n) - double precision, intent(out) :: C(LDC,m) +!END_DOC +!integer, intent(in) :: n,m, LDC, LDP +!double precision, intent(in) :: P(LDP,n) +!double precision, intent(out) :: C(LDC,m) - integer :: info - integer :: i,k - integer :: ipiv(n) - double precision:: tol - double precision, allocatable :: W(:,:), work(:) +!integer :: info +!integer :: i,k +!integer :: ipiv(n) +!double precision:: tol +!double precision, allocatable :: W(:,:), work(:) - allocate(W(LDC,n),work(2*n)) - call svd(P,LDP,C,LDC,W,size(W,1),m,n) +!allocate(W(LDC,n),work(2*n)) +!call svd(P,LDP,C,LDC,W,size(W,1),m,n) - deallocate(W,work) -end +!deallocate(W,work) +!end subroutine svd_mo(n,m,P,LDP,C,LDC) implicit none