LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-07-03 01:45:59 +02:00
Code Issues Releases Wiki Activity

beginning to debug dft

Browse Source
This commit is contained in:
Emmanuel Giner 2017-04-13 14:55:10 +02:00
parent 65a7de3182
commit dab0f90731
6 changed files with 230 additions and 131 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

36
plugins/DFT_Utils/grid_density.irp.f
View File

@ -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

96
plugins/MRPT_Utils/energies_cas.irp.f
View File

@ -294,6 +294,7 @@ 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_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_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)]

5
plugins/MRPT_Utils/mrpt_dress.irp.f
View File

@ -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

181
plugins/MRPT_Utils/mrpt_utils.irp.f
View File

@ -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"

7
plugins/MRPT_Utils/psi_active_prov.irp.f
View File

@ -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

32
src/MO_Basis/cholesky_mo.irp.f
View File

@ -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