LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-11-03 20:54:00 +01:00
Code Issues Releases Wiki Activity

Fixed JM-MRPT2

Browse Source
This commit is contained in:
Anthony Scemama 2017-05-31 19:03:49 +02:00
parent 4fe9c8d844
commit 840fe12c9d
9 changed files with 533 additions and 182 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

0
plugins/MRPT/MRPT_Utils.main.irp.f → plugins/MRPT/jm_mrpt2.irp.f
View File

38
plugins/MRPT/mrpt.irp.f
View File

@ -1,38 +0,0 @@
program MRPT
implicit none
BEGIN_DOC
! TODO
END_DOC
print *, ' _/ '
print *, ' -:\_?, _Jm####La '
print *, 'J"(:" > _]#AZ#Z#UUZ##, '
print *, '_,::./ %(|i%12XmX1*1XL _?, '
print *, ' \..\ _\(vmWQwodY+ia%lnL _",/ ( '
print *, ' .:< ]J=mQD?WXn<uQWmmvd, -.-:=!'
print *, ' "{Z jC]QW|=3Zv)Bi3BmXv3 = _7'
print *, ' ]h[Z6)WQ;)jZs]C;|$BZv+, : ./ '
print *, ' -#sJX%$Wmm#ev]hinW#Xi:` c ; '
print *, ' #X#X23###1}vI$WWmX1>|,)nr" '
print *, ' 4XZ#Xov1v}=)vnXAX1nnv;1n" '
print *, ' ]XX#ZXoovvvivnnnlvvo2*i7 '
print *, ' "23Z#1S2oo2XXSnnnoSo2>v" '
print *, ' miX#L -~`""!!1}oSoe|i7 '
print *, ' 4cn#m, v221=|v[ '
print *, ' ]hI3Zma,;..__wXSe=+vo '
print *, ' ]Zov*XSUXXZXZXSe||vo2 '
print *, ' ]Z#><iiii|i||||==vn2( '
print *, ' ]Z#i<ii||+|=||=:{no2[ '
print *, ' ]ZUsiiiiivi|=||=vo22[ '
print *, ' ]XZvlliiIi|i=|+|vooo '
print *, ' =v1llli||||=|||||lii( '
print *, ' ]iillii||||||||=>=|< '
print *, ' -ziiiii||||||+||==+> '
print *, ' -%|+++||=|=+|=|==/ '
print *, ' -a>====+|====-:- '
print *, ' "~,- -- /- '
print *, ' -. )> '
print *, ' .~ +- '
print *, ' . .... : . '
print *, ' -------~ '
print *, ''
end

193
plugins/MRPT_Utils/density_matrix_based.irp.f Normal file
View File

@ -0,0 +1,193 @@
subroutine contrib_1h2p_dm_based(accu)
implicit none
integer :: i_i,i_r,i_v,i_a,i_b
integer :: i,r,v,a,b
integer :: ispin,jspin
integer :: istate
double precision, intent(out) :: accu(N_states)
double precision :: active_int(n_act_orb,2)
double precision :: delta_e(n_act_orb,2,N_states)
double precision :: get_mo_bielec_integral
accu = 0.d0
!do i_i = 1, 1
do i_i = 1, n_inact_orb
i = list_inact(i_i)
! do i_r = 1, 1
do i_r = 1, n_virt_orb
r = list_virt(i_r)
! do i_v = 1, 1
do i_v = 1, n_virt_orb
v = list_virt(i_v)
do i_a = 1, n_act_orb
a = list_act(i_a)
active_int(i_a,1) = get_mo_bielec_integral(i,a,r,v,mo_integrals_map) ! direct
active_int(i_a,2) = get_mo_bielec_integral(i,a,v,r,mo_integrals_map) ! exchange
do istate = 1, N_states
do jspin=1, 2
delta_e(i_a,jspin,istate) = one_anhil(i_a,jspin,istate) &
- fock_virt_total_spin_trace(r,istate) &
- fock_virt_total_spin_trace(v,istate) &
+ fock_core_inactive_total_spin_trace(i,istate)
delta_e(i_a,jspin,istate) = 1.d0/delta_e(i_a,jspin,istate)
enddo
enddo
enddo
do i_a = 1, n_act_orb
a = list_act(i_a)
do i_b = 1, n_act_orb
! do i_b = i_a, i_a
b = list_act(i_b)
do ispin = 1, 2 ! spin of (i --> r)
do jspin = 1, 2 ! spin of (a --> v)
if(ispin == jspin .and. r.le.v)cycle ! condition not to double count
do istate = 1, N_states
if(ispin == jspin)then
accu(istate) += (active_int(i_a,1) - active_int(i_a,2)) * one_body_dm_mo_spin_index(a,b,istate,ispin) &
* (active_int(i_b,1) - active_int(i_b,2)) &
* delta_e(i_a,jspin,istate)
else
accu(istate) += active_int(i_a,1) * one_body_dm_mo_spin_index(a,b,istate,ispin) * delta_e(i_a,ispin,istate) &
* active_int(i_b,1)
endif
enddo
enddo
enddo
enddo
enddo
enddo
enddo
enddo
end
subroutine contrib_2h1p_dm_based(accu)
implicit none
integer :: i_i,i_j,i_v,i_a,i_b
integer :: i,j,v,a,b
integer :: ispin,jspin
integer :: istate
double precision, intent(out) :: accu(N_states)
double precision :: active_int(n_act_orb,2)
double precision :: delta_e(n_act_orb,2,N_states)
double precision :: get_mo_bielec_integral
accu = 0.d0
do i_i = 1, n_inact_orb
i = list_inact(i_i)
do i_j = 1, n_inact_orb
j = list_inact(i_j)
do i_v = 1, n_virt_orb
v = list_virt(i_v)
do i_a = 1, n_act_orb
a = list_act(i_a)
active_int(i_a,1) = get_mo_bielec_integral(i,j,v,a,mo_integrals_map) ! direct
active_int(i_a,2) = get_mo_bielec_integral(i,j,a,v,mo_integrals_map) ! exchange
do istate = 1, N_states
do jspin=1, 2
! delta_e(i_a,jspin,istate) =
!
delta_e(i_a,jspin,istate) = one_creat(i_a,jspin,istate) - fock_virt_total_spin_trace(v,istate) &
+ fock_core_inactive_total_spin_trace(i,istate) &
+ fock_core_inactive_total_spin_trace(j,istate)
delta_e(i_a,jspin,istate) = 1.d0/delta_e(i_a,jspin,istate)
enddo
enddo
enddo
do i_a = 1, n_act_orb
a = list_act(i_a)
do i_b = 1, n_act_orb
! do i_b = i_a, i_a
b = list_act(i_b)
do ispin = 1, 2 ! spin of (i --> v)
do jspin = 1, 2 ! spin of (j --> a)
if(ispin == jspin .and. i.le.j)cycle ! condition not to double count
do istate = 1, N_states
if(ispin == jspin)then
accu(istate) += (active_int(i_a,1) - active_int(i_a,2)) * one_body_dm_dagger_mo_spin_index(a,b,istate,ispin) &
* (active_int(i_b,1) - active_int(i_b,2)) &
* delta_e(i_a,jspin,istate)
else
accu(istate) += active_int(i_a,1) * one_body_dm_dagger_mo_spin_index(a,b,istate,ispin) * delta_e(i_a,ispin,istate) &
* active_int(i_b,1)
endif
enddo
enddo
enddo
enddo
enddo
enddo
enddo
enddo
end
!subroutine contrib_2p_dm_based(accu)
!implicit none
!integer :: i_r,i_v,i_a,i_b,i_c,i_d
!integer :: r,v,a,b,c,d
!integer :: ispin,jspin
!integer :: istate
!double precision, intent(out) :: accu(N_states)
!double precision :: active_int(n_act_orb,n_act_orb,2)
!double precision :: delta_e(n_act_orb,n_act_orb,2,2,N_states)
!double precision :: get_mo_bielec_integral
!accu = 0.d0
!do i_r = 1, n_virt_orb
! r = list_virt(i_r)
! do i_v = 1, n_virt_orb
! v = list_virt(i_v)
! do i_a = 1, n_act_orb
! a = list_act(i_a)
! do i_b = 1, n_act_orb
! b = list_act(i_b)
! active_int(i_a,i_b,1) = get_mo_bielec_integral(a,b,r,v,mo_integrals_map) ! direct
! active_int(i_a,i_b,2) = get_mo_bielec_integral(a,b,v,r,mo_integrals_map) ! direct
! do istate = 1, N_states
! do jspin=1, 2 ! spin of i_a
! do ispin = 1, 2 ! spin of i_b
! delta_e(i_a,i_b,jspin,ispin,istate) = two_anhil(i_a,i_b,jspin,ispin,istate) &
! - fock_virt_total_spin_trace(r,istate) &
! - fock_virt_total_spin_trace(v,istate)
! delta_e(i_a,i_b,jspin,ispin,istate) = 1.d0/delta_e(i_a,i_b,jspin,ispin,istate)
! enddo
! enddo
! enddo
! enddo
! enddo
! ! diagonal terms
! do i_a = 1, n_act_orb
! a = list_act(i_a)
! do i_b = 1, n_act_orb
! b = list_act(i_b)
! do ispin = 1, 2 ! spin of (a --> r)
! do jspin = 1, 2 ! spin of (b --> v)
! if(ispin == jspin .and. r.le.v)cycle ! condition not to double count
! if(ispin == jspin .and. a.le.b)cycle ! condition not to double count
! do istate = 1, N_states
! if(ispin == jspin)then
! double precision :: contrib_spin
! if(ispin == 1)then
! contrib_spin = two_body_dm_aa_diag_act(i_a,i_b)
! else
! contrib_spin = two_body_dm_bb_diag_act(i_a,i_b)
! endif
! accu(istate) += (active_int(i_a,i_b,1) - active_int(i_a,i_b,2)) * contrib_spin &
! * (active_int(i_a,i_b,1) - active_int(i_a,i_b,2)) &
! * delta_e(i_a,i_b,ispin,jspin,istate)
! else
! accu(istate) += 0.5d0 * active_int(i_a,i_b,1) * two_body_dm_ab_diag_act(i_a,i_b) * delta_e(i_a,i_b,ispin,jspin,istate) &
! * active_int(i_a,i_b,1)
! endif
! enddo
! enddo
! enddo
! enddo
! enddo
! enddo
! enddo
!end

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

@ -1,9 +1,9 @@
BEGIN_PROVIDER [ double precision, energy_cas_dyall, (N_states)]
implicit none
integer :: i
double precision :: energies(N_states_diag)
double precision :: energies(N_states)
do i = 1, N_states
call u0_H_dyall_u0(energies,psi_active,psi_coef,n_det,psi_det_size,psi_det_size,N_states_diag,i)
call u0_H_dyall_u0(energies,psi_active,psi_coef,n_det,psi_det_size,psi_det_size,N_states,i)
energy_cas_dyall(i) = energies(i)
print*, 'energy_cas_dyall(i)', energy_cas_dyall(i)
enddo
@ -13,9 +13,9 @@ END_PROVIDER
BEGIN_PROVIDER [ double precision, energy_cas_dyall_no_exchange, (N_states)]
implicit none
integer :: i
double precision :: energies(N_states_diag)
double precision :: energies(N_states)
do i = 1, N_states
call u0_H_dyall_u0_no_exchange(energies,psi_active,psi_coef,n_det,psi_det_size,psi_det_size,N_states_diag,i)
call u0_H_dyall_u0_no_exchange(energies,psi_active,psi_coef,n_det,psi_det_size,psi_det_size,N_states,i)
energy_cas_dyall_no_exchange(i) = energies(i)
print*, 'energy_cas_dyall(i)_no_exchange', energy_cas_dyall_no_exchange(i)
enddo
@ -28,22 +28,22 @@ BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2,N_states)]
integer :: i,j
integer :: ispin
integer :: orb, hole_particle,spin_exc
double precision :: norm_out(N_states_diag)
double precision :: norm_out(N_states)
integer(bit_kind), allocatable :: psi_in_out(:,:,:)
double precision, allocatable :: psi_in_out_coef(:,:)
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states_diag))
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states))
use bitmasks
integer :: iorb
integer :: state_target
double precision :: energies(n_states_diag)
double precision :: energies(n_states)
do iorb = 1,n_act_orb
do ispin = 1,2
orb = list_act(iorb)
hole_particle = 1
spin_exc = ispin
do i = 1, n_det
do j = 1, n_states_diag
do j = 1, n_states
psi_in_out_coef(i,j) = psi_coef(i,j)
enddo
do j = 1, N_int
@ -53,8 +53,8 @@ BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2,N_states)]
enddo
do state_target = 1,N_states
call apply_exc_to_psi(orb,hole_particle,spin_exc, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target)
one_creat(iorb,ispin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
enddo
enddo
@ -68,22 +68,22 @@ BEGIN_PROVIDER [ double precision, one_anhil, (n_act_orb,2,N_states)]
integer :: i,j
integer :: ispin
integer :: orb, hole_particle,spin_exc
double precision :: norm_out(N_states_diag)
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,N_states_diag))
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states))
integer :: iorb
integer :: state_target
double precision :: energies(n_states_diag)
double precision :: energies(n_states)
do iorb = 1,n_act_orb
do ispin = 1,2
orb = list_act(iorb)
hole_particle = -1
spin_exc = ispin
do i = 1, n_det
do j = 1, n_states_diag
do j = 1, n_states
psi_in_out_coef(i,j) = psi_coef(i,j)
enddo
do j = 1, N_int
@ -93,8 +93,8 @@ BEGIN_PROVIDER [ double precision, one_anhil, (n_act_orb,2,N_states)]
enddo
do state_target = 1, N_states
call apply_exc_to_psi(orb,hole_particle,spin_exc, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target)
one_anhil(iorb,ispin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
enddo
enddo
@ -109,15 +109,15 @@ BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2,N_states)
integer :: ispin,jspin
integer :: orb_i, hole_particle_i,spin_exc_i
integer :: orb_j, hole_particle_j,spin_exc_j
double precision :: norm_out(N_states_diag)
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,N_states_diag))
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states))
integer :: iorb,jorb
integer :: state_target
double precision :: energies(n_states_diag)
double precision :: energies(n_states)
do iorb = 1,n_act_orb
do ispin = 1,2
orb_i = list_act(iorb)
@ -129,7 +129,7 @@ BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2,N_states)
hole_particle_j = 1
spin_exc_j = jspin
do i = 1, n_det
do j = 1, n_states_diag
do j = 1, n_states
psi_in_out_coef(i,j) = psi_coef(i,j)
enddo
do j = 1, N_int
@ -139,10 +139,10 @@ BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2,N_states)
enddo
do state_target = 1 , N_states
call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target)
two_creat(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
enddo
enddo
@ -159,16 +159,16 @@ BEGIN_PROVIDER [ double precision, two_anhil, (n_act_orb,n_act_orb,2,2,N_states)
integer :: ispin,jspin
integer :: orb_i, hole_particle_i,spin_exc_i
integer :: orb_j, hole_particle_j,spin_exc_j
double precision :: norm_out(N_states_diag)
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,N_states_diag))
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states))
integer :: iorb,jorb
integer :: state_target
state_target = 1
double precision :: energies(n_states_diag)
double precision :: energies(n_states)
do iorb = 1,n_act_orb
do ispin = 1,2
orb_i = list_act(iorb)
@ -180,7 +180,7 @@ BEGIN_PROVIDER [ double precision, two_anhil, (n_act_orb,n_act_orb,2,2,N_states)
hole_particle_j = -1
spin_exc_j = jspin
do i = 1, n_det
do j = 1, n_states_diag
do j = 1, n_states
psi_in_out_coef(i,j) = psi_coef(i,j)
enddo
do j = 1, N_int
@ -189,10 +189,10 @@ BEGIN_PROVIDER [ double precision, two_anhil, (n_act_orb,n_act_orb,2,2,N_states)
enddo
enddo
call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target)
two_anhil(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
enddo
enddo
@ -208,15 +208,15 @@ BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2
integer :: ispin,jspin
integer :: orb_i, hole_particle_i,spin_exc_i
integer :: orb_j, hole_particle_j,spin_exc_j
double precision :: norm_out(N_states_diag)
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,N_states_diag))
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states))
integer :: iorb,jorb
integer :: state_target
double precision :: energies(n_states_diag)
double precision :: energies(n_states)
do iorb = 1,n_act_orb
do ispin = 1,2
orb_i = list_act(iorb)
@ -228,7 +228,7 @@ BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2
hole_particle_j = -1
spin_exc_j = jspin
do i = 1, n_det
do j = 1, n_states_diag
do j = 1, n_states
psi_in_out_coef(i,j) = psi_coef(i,j)
enddo
do j = 1, N_int
@ -238,14 +238,14 @@ BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2
enddo
do state_target = 1, N_states
call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
if(orb_i == orb_j .and. ispin .ne. jspin)then
call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target)
one_anhil_one_creat(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall_no_exchange(state_target) - energies(state_target)
else
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target)
one_anhil_one_creat(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
endif
enddo
@ -264,16 +264,16 @@ BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_a
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_diag)
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,N_states_diag))
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states))
integer :: iorb,jorb
integer :: korb
integer :: state_target
double precision :: energies(n_states_diag)
double precision :: energies(n_states)
do iorb = 1,n_act_orb
do ispin = 1,2
orb_i = list_act(iorb)
@ -290,7 +290,7 @@ BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_a
hole_particle_k = -1
spin_exc_k = kspin
do i = 1, n_det
do j = 1, n_states_diag
do j = 1, n_states
psi_in_out_coef(i,j) = psi_coef(i,j)
enddo
do j = 1, N_int
@ -301,12 +301,12 @@ BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_a
do state_target = 1, N_states
call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
call apply_exc_to_psi(orb_k,hole_particle_k,spin_exc_k, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target)
two_anhil_one_creat(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
enddo
enddo
@ -326,16 +326,16 @@ BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_a
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_diag)
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,N_states_diag))
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states))
integer :: iorb,jorb
integer :: korb
integer :: state_target
double precision :: energies(n_states_diag)
double precision :: energies(n_states)
do iorb = 1,n_act_orb
do ispin = 1,2
orb_i = list_act(iorb)
@ -352,7 +352,7 @@ BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_a
hole_particle_k = -1
spin_exc_k = kspin
do i = 1, n_det
do j = 1, n_states_diag
do j = 1, n_states
psi_in_out_coef(i,j) = psi_coef(i,j)
enddo
do j = 1, N_int
@ -362,12 +362,12 @@ BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_a
enddo
do state_target = 1, N_states
call apply_exc_to_psi(orb_k,hole_particle_k,spin_exc_k, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target)
two_creat_one_anhil(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
enddo
enddo
@ -387,16 +387,16 @@ BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2
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_diag)
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,N_states_diag))
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states))
integer :: iorb,jorb
integer :: korb
integer :: state_target
double precision :: energies(n_states_diag)
double precision :: energies(n_states)
do iorb = 1,n_act_orb
do ispin = 1,2
orb_i = list_act(iorb)
@ -413,7 +413,7 @@ BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2
hole_particle_k = 1
spin_exc_k = kspin
do i = 1, n_det
do j = 1, n_states_diag
do j = 1, n_states
psi_in_out_coef(i,j) = psi_coef(i,j)
enddo
do j = 1, N_int
@ -423,12 +423,12 @@ BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2
enddo
do state_target = 1, N_states
call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
call apply_exc_to_psi(orb_k,hole_particle_k,spin_exc_k, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target)
three_creat(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
enddo
enddo
@ -448,16 +448,16 @@ BEGIN_PROVIDER [ double precision, three_anhil, (n_act_orb,n_act_orb,n_act_orb,2
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_diag)
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,N_states_diag))
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states))
integer :: iorb,jorb
integer :: korb
integer :: state_target
double precision :: energies(n_states_diag)
double precision :: energies(n_states)
do iorb = 1,n_act_orb
do ispin = 1,2
orb_i = list_act(iorb)
@ -474,7 +474,7 @@ BEGIN_PROVIDER [ double precision, three_anhil, (n_act_orb,n_act_orb,n_act_orb,2
hole_particle_k = -1
spin_exc_k = kspin
do i = 1, n_det
do j = 1, n_states_diag
do j = 1, n_states
psi_in_out_coef(i,j) = psi_coef(i,j)
enddo
do j = 1, N_int
@ -484,12 +484,12 @@ BEGIN_PROVIDER [ double precision, three_anhil, (n_act_orb,n_act_orb,n_act_orb,2
enddo
do state_target = 1, N_states
call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
call apply_exc_to_psi(orb_k,hole_particle_k,spin_exc_k, &
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target)
three_anhil(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
enddo
enddo
@ -511,15 +511,15 @@ END_PROVIDER
integer :: ispin,jspin
integer :: orb_i, hole_particle_i
integer :: orb_v
double precision :: norm_out(N_states_diag)
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,N_states_diag))
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states))
integer :: iorb,jorb,i_ok
integer :: state_target
double precision :: energies(n_states_diag)
double precision :: energies(n_states)
double precision :: hij
double precision :: norm(N_states,2),norm_no_inv(N_states,2),norm_bis(N_states,2)
double precision :: energies_alpha_beta(N_states,2)
@ -585,7 +585,7 @@ END_PROVIDER
energies_alpha_beta(state_target, ispin) = - mo_bielec_integral_jj_exchange(orb_i,orb_v)
! energies_alpha_beta(state_target, ispin) = 0.d0
if(norm(state_target,ispin) .ne. 0.d0 .and. dabs(norm_no_inv(state_target,ispin)) .gt. thresh_norm)then
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target)
energies_alpha_beta(state_target, ispin) += energies(state_target)
endif
enddo
@ -616,15 +616,15 @@ BEGIN_PROVIDER [ double precision, one_anhil_inact, (n_inact_orb,n_act_orb,N_Sta
integer :: i,iorb,j
integer :: ispin,jspin
integer :: orb_i, hole_particle_i
double precision :: norm_out(N_states_diag)
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,N_states_diag))
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states))
integer :: jorb,i_ok,aorb,orb_a
integer :: state_target
double precision :: energies(n_states_diag)
double precision :: energies(n_states)
double precision :: hij
double precision :: norm(N_states,2),norm_no_inv(N_states,2)
double precision :: energies_alpha_beta(N_states,2)
@ -688,7 +688,7 @@ BEGIN_PROVIDER [ double precision, one_anhil_inact, (n_inact_orb,n_act_orb,N_Sta
do state_target = 1, N_states
energies_alpha_beta(state_target, ispin) = 0.d0
if(norm(state_target,ispin) .ne. 0.d0 .and. dabs(norm_no_inv(state_target,ispin)) .gt. thresh_norm)then
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target)
energies_alpha_beta(state_target, ispin) += energies(state_target)
endif
enddo
@ -701,11 +701,6 @@ BEGIN_PROVIDER [ double precision, one_anhil_inact, (n_inact_orb,n_act_orb,N_Sta
else
one_anhil_inact(iorb,aorb,state_target) = 0.d0
endif
! print*, '********'
! print*, energies_alpha_beta(state_target,1) , energies_alpha_beta(state_target,2)
! print*, norm_bis(state_target,1) , norm_bis(state_target,2)
! print*, one_anhil_inact(iorb,aorb,state_target)
! print*, one_creat(aorb,1,state_target)
enddo
enddo
enddo
@ -719,15 +714,15 @@ BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_State
integer :: ispin,jspin
integer :: orb_i, hole_particle_i
integer :: orb_v
double precision :: norm_out(N_states_diag)
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,N_states_diag))
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states))
integer :: iorb,jorb,i_ok,aorb,orb_a
integer :: state_target
double precision :: energies(n_states_diag)
double precision :: energies(n_states)
double precision :: hij
double precision :: norm(N_states,2),norm_no_inv(N_states,2)
double precision :: energies_alpha_beta(N_states,2)
@ -791,7 +786,7 @@ BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_State
do state_target = 1, N_states
energies_alpha_beta(state_target, ispin) = 0.d0
if(norm(state_target,ispin) .ne. 0.d0 .and. dabs(norm_no_inv(state_target,ispin)) .gt. thresh_norm)then
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target)
! print*, energies(state_target)
energies_alpha_beta(state_target, ispin) += energies(state_target)
endif
@ -825,19 +820,19 @@ END_PROVIDER
integer :: ispin,jspin
integer :: orb_i, hole_particle_i
integer :: orb_v
double precision :: norm_out(N_states_diag),diag_elem(N_det),interact_psi0(N_det)
double precision :: norm_out(N_states),diag_elem(N_det),interact_psi0(N_det)
double precision :: delta_e_inact_virt(N_states)
integer(bit_kind), allocatable :: psi_in_out(:,:,:)
double precision, allocatable :: psi_in_out_coef(:,:)
double precision, allocatable :: H_matrix(:,:),eigenvectors(:,:),eigenvalues(:)
use bitmasks
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states_diag),H_matrix(N_det+1,N_det+1))
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states),H_matrix(N_det+1,N_det+1))
allocate (eigenvectors(size(H_matrix,1),N_det+1))
allocate (eigenvalues(N_det+1))
integer :: iorb,jorb,i_ok
integer :: state_target
double precision :: energies(n_states_diag)
double precision :: energies(n_states)
double precision :: hij
double precision :: energies_alpha_beta(N_states,2)
@ -973,21 +968,21 @@ subroutine give_singles_and_partial_doubles_1h1p_contrib(matrix_1h1p,e_corr_from
integer :: ispin,jspin
integer :: orb_i, hole_particle_i
integer :: orb_v
double precision :: norm_out(N_states_diag),diag_elem(N_det),interact_psi0(N_det)
double precision :: norm_out(N_states),diag_elem(N_det),interact_psi0(N_det)
double precision :: delta_e_inact_virt(N_states)
integer(bit_kind), allocatable :: psi_in_out(:,:,:)
double precision, allocatable :: psi_in_out_coef(:,:)
double precision, allocatable :: H_matrix(:,:),eigenvectors(:,:),eigenvalues(:),interact_cas(:,:)
double precision, allocatable :: delta_e_det(:,:)
use bitmasks
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states_diag),H_matrix(N_det+1,N_det+1))
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states),H_matrix(N_det+1,N_det+1))
allocate (eigenvectors(size(H_matrix,1),N_det+1))
allocate (eigenvalues(N_det+1),interact_cas(N_det,N_det))
allocate (delta_e_det(N_det,N_det))
integer :: iorb,jorb,i_ok
integer :: state_target
double precision :: energies(n_states_diag)
double precision :: energies(n_states)
double precision :: hij
double precision :: energies_alpha_beta(N_states,2)
double precision :: lamda_pt2(N_det)
@ -1023,7 +1018,7 @@ subroutine give_singles_and_partial_doubles_1h1p_contrib(matrix_1h1p,e_corr_from
interact_psi0(i) = 0.d0
do j = 1 , N_det
call i_H_j(psi_in_out(1,1,i),psi_det(1,1,j),N_int,hij)
call get_delta_e_dyall(psi_det(1,1,j),psi_in_out(1,1,i),coef_array,hij,delta_e_det(i,j))
call get_delta_e_dyall(psi_det(1,1,j),psi_in_out(1,1,i),delta_e_det(i,j))
interact_cas(i,j) = hij
interact_psi0(i) += hij * psi_coef(j,1)
enddo

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

@ -84,7 +84,8 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip
do i_state = 1, N_states
coef_array(i_state) = psi_coef(index_i,i_state)
enddo
call get_delta_e_dyall(psi_det(1,1,index_i),tq(1,1,i_alpha),coef_array,hialpha,delta_e)
call get_delta_e_dyall(psi_det(1,1,index_i),tq(1,1,i_alpha),delta_e)
! call get_delta_e_dyall_general_mp(psi_det(1,1,index_i),tq(1,1,i_alpha),delta_e)
hij_array(index_i) = hialpha
call get_excitation(psi_det(1,1,index_i),tq(1,1,i_alpha),exc,degree,phase,N_int)
! phase_array(index_i) = phase

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

@ -58,8 +58,6 @@
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
@ -121,7 +119,7 @@
! 1h2p
delta_ij_tmp = 0.d0
!call give_1h2p_contrib(delta_ij_tmp)
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
@ -137,7 +135,7 @@
! 2h1p
delta_ij_tmp = 0.d0
!call give_2h1p_contrib(delta_ij_tmp)
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
@ -223,9 +221,9 @@ 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_eigenvectors, (N_det,N_states_diag) ]
&BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_eigenvectors_s2, (N_states_diag) ]
BEGIN_PROVIDER [ double precision, CI_electronic_dressed_pt2_new_energy, (N_states) ]
&BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_eigenvectors, (N_det,N_states) ]
&BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_eigenvectors_s2, (N_states) ]
BEGIN_DOC
! Eigenvectors/values of the CI matrix
END_DOC
@ -244,14 +242,14 @@ END_PROVIDER
double precision, allocatable :: e_array(:)
integer, allocatable :: iorder(:)
! Guess values for the "N_states_diag" states of the CI_dressed_pt2_new_eigenvectors
do j=1,min(N_states_diag,N_det)
! Guess values for the "N_states" states of the CI_dressed_pt2_new_eigenvectors
do j=1,min(N_states,N_det)
do i=1,N_det
CI_dressed_pt2_new_eigenvectors(i,j) = psi_coef(i,j)
enddo
enddo
do j=N_det+1,N_states_diag
do j=N_det+1,N_states
do i=1,N_det
CI_dressed_pt2_new_eigenvectors(i,j) = 0.d0
enddo
@ -267,8 +265,8 @@ 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)
CI_electronic_energy(:) = 0.d0
Hmatrix_dressed_pt2_new_symmetrized(1,1,1),size(H_matrix_all_dets,1),N_det)
CI_electronic_dressed_pt2_new_energy(:) = 0.d0
if (s2_eig) then
i_state = 0
allocate (s2_eigvalues(N_det))
@ -291,54 +289,54 @@ 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_electronic_dressed_pt2_new_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
if(good_state_array(j))cycle
i_other_state +=1
if(i_state+i_other_state.gt.n_states_diag)then
if(i_state+i_other_state.gt.n_states)then
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_electronic_dressed_pt2_new_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
print*,''
print*,'!!!!!!!! WARNING !!!!!!!!!'
print*,' Within the ',N_det,'determinants selected'
print*,' and the ',N_states_diag,'states requested'
print*,' and the ',N_states,'states requested'
print*,' We did not find any state with S^2 values close to ',expected_s2
print*,' We will then set the first N_states eigenvectors of the H matrix'
print*,' as the CI_eigenvectors'
print*,' as the CI_dressed_pt2_new_eigenvectors'
print*,' You should consider more states and maybe ask for s2_eig to be .True. or just enlarge the CI space'
print*,''
do j=1,min(N_states_diag,N_det)
do j=1,min(N_states,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_electronic_dressed_pt2_new_energy(j) = eigenvalues(j)
CI_dressed_pt2_new_eigenvectors_s2(j) = s2_eigvalues(j)
enddo
endif
deallocate(index_good_state_array,good_state_array)
deallocate(s2_eigvalues)
else
call u_0_S2_u_0(CI_eigenvectors_s2,eigenvectors,N_det,psi_det,N_int,&
min(N_det,N_states_diag),size(eigenvectors,1))
! Select the "N_states_diag" states of lowest energy
do j=1,min(N_det,N_states_diag)
call u_0_S2_u_0(CI_dressed_pt2_new_eigenvectors_s2,eigenvectors,N_det,psi_det,N_int,&
min(N_det,N_states),size(eigenvectors,1))
! Select the "N_states" states of lowest energy
do j=1,min(N_det,N_states)
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_electronic_dressed_pt2_new_energy(j) = eigenvalues(j)
enddo
endif
deallocate(eigenvectors,eigenvalues)
@ -348,7 +346,7 @@ END_PROVIDER
END_PROVIDER
BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_energy, (N_states_diag) ]
BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_energy, (N_states) ]
implicit none
BEGIN_DOC
! N_states lowest eigenvalues of the CI matrix
@ -357,11 +355,11 @@ BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_energy, (N_states_diag) ]
integer :: j
character*(8) :: st
call write_time(output_determinants)
do j=1,N_states_diag
do j=1,N_states
CI_dressed_pt2_new_energy(j) = CI_electronic_dressed_pt2_new_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))
call write_double(output_determinants,CI_dressed_pt2_new_eigenvectors_s2(j),'S^2 of state '//trim(st))
enddo
END_PROVIDER

2
plugins/MRPT_Utils/new_way.irp.f
View File

@ -804,7 +804,7 @@ subroutine give_1p_sec_order_singles_contrib(matrix_1p)
do jdet = 1,N_det
double precision :: coef_array(N_states),hij_test
call i_H_j(det_tmp,psi_det(1,1,jdet),N_int,himono)
call get_delta_e_dyall(psi_det(1,1,jdet),det_tmp,coef_array,hij_test,delta_e)
call get_delta_e_dyall(psi_det(1,1,jdet),det_tmp,delta_e)
do state_target = 1, N_states
! matrix_1p(idet,jdet,state_target) += himono * coef_det_pert(i,r,ispin,state_target,1)
matrix_1p(idet,jdet,state_target) += himono * hij_det_pert(i,r,ispin) / delta_e(state_target)

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

@ -152,7 +152,7 @@ subroutine give_particles_in_virt_space(det_1,n_particles_spin,n_particles,parti
end
subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final)
subroutine get_delta_e_dyall(det_1,det_2,delta_e_final)
BEGIN_DOC
! routine that returns the delta_e with the Moller Plesset and Dyall operators
!
@ -170,7 +170,6 @@ subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final)
use bitmasks
double precision, intent(out) :: delta_e_final(N_states)
integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2)
double precision, intent(in) :: coef_array(N_states),hij
integer :: i,j,k,l
integer :: i_state
@ -293,20 +292,9 @@ subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final)
if (n_holes_act == 0 .and. n_particles_act == 1) then
ispin = particle_list_practical(1,1)
i_particle_act = particle_list_practical(2,1)
! call get_excitation_degree(det_1,det_2,degree,N_int)
! if(degree == 1)then
! call get_excitation(det_1,det_2,exc,degree,phase,N_int)
! call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
! i_hole = list_inact_reverse(h1)
! i_part = list_act_reverse(p1)
! do i_state = 1, N_states
! delta_e_act(i_state) += one_anhil_inact(i_hole,i_part,i_state)
! enddo
! else if (degree == 2)then
do i_state = 1, N_states
delta_e_act(i_state) += one_creat(i_particle_act,ispin,i_state)
enddo
! endif
else if (n_holes_act == 1 .and. n_particles_act == 0) then
ispin = hole_list_practical(1,1)
@ -433,3 +421,159 @@ subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final)
end
subroutine get_delta_e_dyall_general_mp(det_1,det_2,delta_e_final)
BEGIN_DOC
! routine that returns the delta_e with the Moller Plesset and Dyall operators
!
! with det_1 being a determinant from the cas, and det_2 being a perturber
!
! Delta_e(det_1,det_2) = sum (hole) epsilon(hole) + sum(part) espilon(part) + delta_e(act)
!
! where hole is necessary in the inactive, part necessary in the virtuals
!
! and delta_e(act) is obtained as the sum of energies of excitations a la MP
!
END_DOC
implicit none
use bitmasks
double precision, intent(out) :: delta_e_final(N_states)
integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2)
integer :: i,j,k,l
integer :: i_state
integer :: n_holes_spin(2)
integer :: n_holes
integer :: holes_list(N_int*bit_kind_size,2)
double precision :: delta_e_inactive(N_states)
integer :: i_hole_inact
call give_holes_in_inactive_space(det_2,n_holes_spin,n_holes,holes_list)
delta_e_inactive = 0.d0
do i = 1, n_holes_spin(1)
i_hole_inact = holes_list(i,1)
do i_state = 1, N_states
delta_e_inactive += fock_core_inactive_total_spin_trace(i_hole_inact,i_state)
enddo
enddo
do i = 1, n_holes_spin(2)
i_hole_inact = holes_list(i,2)
do i_state = 1, N_states
delta_e_inactive(i_state) += fock_core_inactive_total_spin_trace(i_hole_inact,i_state)
enddo
enddo
double precision :: delta_e_virt(N_states)
integer :: i_part_virt
integer :: n_particles_spin(2)
integer :: n_particles
integer :: particles_list(N_int*bit_kind_size,2)
call give_particles_in_virt_space(det_2,n_particles_spin,n_particles,particles_list)
delta_e_virt = 0.d0
do i = 1, n_particles_spin(1)
i_part_virt = particles_list(i,1)
do i_state = 1, N_states
delta_e_virt += fock_virt_total_spin_trace(i_part_virt,i_state)
enddo
enddo
do i = 1, n_particles_spin(2)
i_part_virt = particles_list(i,2)
do i_state = 1, N_states
delta_e_virt += fock_virt_total_spin_trace(i_part_virt,i_state)
enddo
enddo
integer :: n_holes_spin_act(2),n_particles_spin_act(2)
integer :: n_holes_act,n_particles_act
integer :: holes_active_list(2*n_act_orb,2)
integer :: holes_active_list_spin_traced(4*n_act_orb)
integer :: particles_active_list(2*n_act_orb,2)
integer :: particles_active_list_spin_traced(4*n_act_orb)
double precision :: delta_e_act(N_states)
delta_e_act = 0.d0
call give_holes_and_particles_in_active_space(det_1,det_2,n_holes_spin_act,n_particles_spin_act, &
n_holes_act,n_particles_act,holes_active_list,particles_active_list)
integer :: icount,icountbis
integer :: hole_list_practical(2,elec_num_tab(1)+elec_num_tab(2)), particle_list_practical(2,elec_num_tab(1)+elec_num_tab(2))
icount = 0
icountbis = 0
do i = 1, n_holes_spin_act(1)
icount += 1
icountbis += 1
hole_list_practical(1,icountbis) = 1 ! spin
hole_list_practical(2,icountbis) = holes_active_list(i,1) ! index of active orb
holes_active_list_spin_traced(icount) = holes_active_list(i,1)
enddo
do i = 1, n_holes_spin_act(2)
icount += 1
icountbis += 1
hole_list_practical(1,icountbis) = 2
hole_list_practical(2,icountbis) = holes_active_list(i,2)
holes_active_list_spin_traced(icount) = holes_active_list(i,2)
enddo
if(icount .ne. n_holes_act) then
print*,''
print*, icount, n_holes_act
print * , 'pb in holes_active_list_spin_traced !!'
stop
endif
icount = 0
icountbis = 0
do i = 1, n_particles_spin_act(1)
icount += 1
icountbis += 1
particle_list_practical(1,icountbis) = 1
particle_list_practical(2,icountbis) = particles_active_list(i,1)
particles_active_list_spin_traced(icount) = particles_active_list(i,1)
enddo
do i = 1, n_particles_spin_act(2)
icount += 1
icountbis += 1
particle_list_practical(1,icountbis) = 2
particle_list_practical(2,icountbis) = particles_active_list(i,2)
particles_active_list_spin_traced(icount) = particles_active_list(i,2)
enddo
if(icount .ne. n_particles_act) then
print*, icount, n_particles_act
print * , 'pb in particles_active_list_spin_traced !!'
stop
endif
integer :: i_hole_act, j_hole_act, k_hole_act
integer :: i_particle_act, j_particle_act, k_particle_act
integer :: ispin,jspin,kspin
do i = 1, n_holes_act
ispin = hole_list_practical(1,i)
i_hole_act = hole_list_practical(2,i)
do i_state = 1, N_states
delta_e_act(i_state) += one_anhil(i_hole_act , ispin,i_state)
enddo
enddo
do i = 1, n_particles_act
ispin = particle_list_practical(1,i)
i_particle_act = particle_list_practical(2,i)
do i_state = 1, N_states
delta_e_act(i_state) += one_creat(i_particle_act, ispin,i_state)
enddo
enddo
do i_state = 1, n_states
delta_e_final(i_state) = delta_e_act(i_state) + delta_e_inactive(i_state) - delta_e_virt(i_state)
enddo
end

60
src/Determinants/density_matrix.irp.f
View File

@ -15,6 +15,57 @@
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, one_body_dm_mo_spin_index, (mo_tot_num_align,mo_tot_num,N_states,2) ]
implicit none
integer :: i,j,ispin,istate
ispin = 1
do istate = 1, N_states
do j = 1, mo_tot_num
do i = 1, mo_tot_num
one_body_dm_mo_spin_index(i,j,istate,ispin) = one_body_dm_mo_alpha(i,j,istate)
enddo
enddo
enddo
ispin = 2
do istate = 1, N_states
do j = 1, mo_tot_num
do i = 1, mo_tot_num
one_body_dm_mo_spin_index(i,j,istate,ispin) = one_body_dm_mo_beta(i,j,istate)
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, one_body_dm_dagger_mo_spin_index, (mo_tot_num_align,mo_tot_num,N_states,2) ]
implicit none
integer :: i,j,ispin,istate
ispin = 1
do istate = 1, N_states
do j = 1, mo_tot_num
one_body_dm_dagger_mo_spin_index(j,j,istate,ispin) = 1 - one_body_dm_mo_alpha(j,j,istate)
do i = j+1, mo_tot_num
one_body_dm_dagger_mo_spin_index(i,j,istate,ispin) = -one_body_dm_mo_alpha(i,j,istate)
one_body_dm_dagger_mo_spin_index(j,i,istate,ispin) = -one_body_dm_mo_alpha(i,j,istate)
enddo
enddo
enddo
ispin = 2
do istate = 1, N_states
do j = 1, mo_tot_num
one_body_dm_dagger_mo_spin_index(j,j,istate,ispin) = 1 - one_body_dm_mo_beta(j,j,istate)
do i = j+1, mo_tot_num
one_body_dm_dagger_mo_spin_index(i,j,istate,ispin) = -one_body_dm_mo_beta(i,j,istate)
one_body_dm_dagger_mo_spin_index(j,i,istate,ispin) = -one_body_dm_mo_beta(i,j,istate)
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, one_body_dm_mo_alpha, (mo_tot_num_align,mo_tot_num,N_states) ]
&BEGIN_PROVIDER [ double precision, one_body_dm_mo_beta, (mo_tot_num_align,mo_tot_num,N_states) ]
implicit none
@ -285,6 +336,8 @@ END_PROVIDER
BEGIN_PROVIDER [ double precision, one_body_dm_ao_alpha, (ao_num_align,ao_num) ]
&BEGIN_PROVIDER [ double precision, one_body_dm_ao_beta, (ao_num_align,ao_num) ]
&BEGIN_PROVIDER [ double precision, one_body_dm_ao_alpha_no_align, (ao_num,ao_num) ]
&BEGIN_PROVIDER [ double precision, one_body_dm_ao_beta_no_align, (ao_num,ao_num) ]
BEGIN_DOC
! one body density matrix on the AO basis : rho_AO(alpha) , rho_AO(beta)
END_DOC
@ -303,11 +356,16 @@ END_PROVIDER
! if(dabs(dm_mo).le.1.d-10)cycle
one_body_dm_ao_alpha(l,k) += mo_coef(k,i) * mo_coef(l,j) * mo_alpha
one_body_dm_ao_beta(l,k) += mo_coef(k,i) * mo_coef(l,j) * mo_beta
enddo
enddo
enddo
enddo
do i = 1, ao_num
do j = 1, ao_num
one_body_dm_ao_alpha_no_align(j,i) = one_body_dm_ao_alpha(j,i)
one_body_dm_ao_beta_no_align(j,i) = one_body_dm_ao_beta(j,i)
enddo
enddo
END_PROVIDER