mirror of
https://github.com/LCPQ/quantum_package
synced 2024-12-23 04:43:50 +01:00
mrpt new with multi state version
This commit is contained in:
Emmanuel Giner
parent
6ebeae0a10
commit
dbf894a99a
@ -18,7 +18,7 @@ IRPF90_FLAGS : --ninja --align=32
|
||||
# 0 : Deactivate
|
||||
#
|
||||
[OPTION]
|
||||
MODE : OPT ; [ OPT | PROFILE | DEBUG ] : Chooses the section below
|
||||
MODE : DEBUG ; [ OPT | PROFILE | DEBUG ] : Chooses the section below
|
||||
CACHE : 1 ; Enable cache_compile.py
|
||||
OPENMP : 1 ; Append OpenMP flags
|
||||
|
||||
|
||||
@ -27,43 +27,17 @@ subroutine routine_2
|
||||
implicit none
|
||||
integer :: i
|
||||
do i = 1, n_core_inact_orb
|
||||
print*,fock_core_inactive_total(i,1),fock_core_inactive(i)
|
||||
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)
|
||||
print*,one_body_dm_mo_alpha(j_act_orb,j_act_orb),one_body_dm_mo_beta(j_act_orb,j_act_orb)
|
||||
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
|
||||
|
||||
end
|
||||
|
||||
subroutine routine
|
||||
implicit none
|
||||
integer :: i,j
|
||||
integer :: orb, spin_exc
|
||||
integer :: hole_particle
|
||||
double precision, allocatable :: norm_out(:)
|
||||
allocate(norm_out(N_states_diag))
|
||||
|
||||
orb = list_virt(10)
|
||||
hole_particle = -1
|
||||
spin_exc = 1
|
||||
|
||||
call apply_exc_to_psi(orb,hole_particle,spin_exc, &
|
||||
norm_out,psi_det,psi_coef, n_det,psi_det_size,psi_det_size,N_states_diag)
|
||||
do i = 1, N_det
|
||||
if(psi_coef(i,1).ne.0.d0)then
|
||||
print*, ''
|
||||
call debug_det(psi_det(1,1,i),N_int)
|
||||
print*, 'coef = ',psi_coef(i,1)
|
||||
endif
|
||||
enddo
|
||||
print*,'norm_out = ',norm_out
|
||||
|
||||
deallocate(norm_out)
|
||||
|
||||
end
|
||||
|
||||
@ -10,7 +10,7 @@ BEGIN_PROVIDER [ double precision, energy_cas_dyall, (N_states)]
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2)]
|
||||
BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2,N_states)]
|
||||
implicit none
|
||||
integer :: i,j
|
||||
integer :: ispin
|
||||
@ -21,6 +21,8 @@ BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2)]
|
||||
use bitmasks
|
||||
|
||||
integer :: iorb
|
||||
integer :: state_target
|
||||
double precision :: energies(n_states_diag)
|
||||
do iorb = 1,n_act_orb
|
||||
do ispin = 1,2
|
||||
orb = list_act(iorb)
|
||||
@ -35,19 +37,18 @@ BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2)]
|
||||
psi_in_out(j,2,i) = psi_det(j,2,i)
|
||||
enddo
|
||||
enddo
|
||||
integer :: state_target
|
||||
state_target = 1
|
||||
double precision :: energies(n_states_diag)
|
||||
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)
|
||||
one_creat(iorb,ispin) = energy_cas_dyall(state_target) - energies(state_target)
|
||||
one_creat(iorb,ispin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, one_anhil, (n_act_orb,2)]
|
||||
BEGIN_PROVIDER [ double precision, one_anhil, (n_act_orb,2,N_states)]
|
||||
implicit none
|
||||
integer :: i,j
|
||||
integer :: ispin
|
||||
@ -58,6 +59,8 @@ BEGIN_PROVIDER [ double precision, one_anhil, (n_act_orb,2)]
|
||||
use bitmasks
|
||||
|
||||
integer :: iorb
|
||||
integer :: state_target
|
||||
double precision :: energies(n_states_diag)
|
||||
do iorb = 1,n_act_orb
|
||||
do ispin = 1,2
|
||||
orb = list_act(iorb)
|
||||
@ -72,30 +75,18 @@ BEGIN_PROVIDER [ double precision, one_anhil, (n_act_orb,2)]
|
||||
psi_in_out(j,2,i) = psi_det(j,2,i)
|
||||
enddo
|
||||
enddo
|
||||
integer :: state_target
|
||||
state_target = 1
|
||||
double precision :: energies(n_states_diag)
|
||||
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)
|
||||
! do j = 1, n_det
|
||||
! print*, 'psi_in_out_coef'
|
||||
! print*, psi_in_out_coef(j,1)
|
||||
! call debug_det(psi_in_out(1,1,j),N_int)
|
||||
! enddo
|
||||
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
|
||||
! print*,'energy_cas_dyall(state_target)'
|
||||
! print*, energy_cas_dyall(state_target)
|
||||
! print*,'energies(state_target)'
|
||||
! print*, energies(state_target)
|
||||
one_anhil(iorb,ispin) = energy_cas_dyall(state_target) - energies(state_target)
|
||||
! print*,'one_anhil(iorb,ispin)'
|
||||
! print*, one_anhil(iorb,ispin)
|
||||
one_anhil(iorb,ispin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2)]
|
||||
BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2,N_states)]
|
||||
implicit none
|
||||
integer :: i,j
|
||||
integer :: ispin,jspin
|
||||
@ -108,7 +99,6 @@ BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2)]
|
||||
|
||||
integer :: iorb,jorb
|
||||
integer :: state_target
|
||||
state_target = 1
|
||||
double precision :: energies(n_states_diag)
|
||||
do iorb = 1,n_act_orb
|
||||
do ispin = 1,2
|
||||
@ -129,12 +119,14 @@ BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2)]
|
||||
psi_in_out(j,2,i) = psi_det(j,2,i)
|
||||
enddo
|
||||
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)
|
||||
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)
|
||||
two_creat(iorb,jorb,ispin,jspin) = energy_cas_dyall(state_target) - energies(state_target)
|
||||
two_creat(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
@ -142,7 +134,7 @@ BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2)]
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, two_anhil, (n_act_orb,n_act_orb,2,2)]
|
||||
BEGIN_PROVIDER [ double precision, two_anhil, (n_act_orb,n_act_orb,2,2,N_states)]
|
||||
implicit none
|
||||
integer :: i,j
|
||||
integer :: ispin,jspin
|
||||
@ -181,7 +173,7 @@ BEGIN_PROVIDER [ double precision, two_anhil, (n_act_orb,n_act_orb,2,2)]
|
||||
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)
|
||||
two_anhil(iorb,jorb,ispin,jspin) = energy_cas_dyall(state_target) - energies(state_target)
|
||||
two_anhil(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
@ -189,7 +181,7 @@ BEGIN_PROVIDER [ double precision, two_anhil, (n_act_orb,n_act_orb,2,2)]
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2)]
|
||||
BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2,N_States)]
|
||||
implicit none
|
||||
integer :: i,j
|
||||
integer :: ispin,jspin
|
||||
@ -202,7 +194,6 @@ BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2
|
||||
|
||||
integer :: iorb,jorb
|
||||
integer :: state_target
|
||||
state_target = 1
|
||||
double precision :: energies(n_states_diag)
|
||||
do iorb = 1,n_act_orb
|
||||
do ispin = 1,2
|
||||
@ -223,12 +214,14 @@ BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2
|
||||
psi_in_out(j,2,i) = psi_det(j,2,i)
|
||||
enddo
|
||||
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)
|
||||
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)
|
||||
one_anhil_one_creat(iorb,jorb,ispin,jspin) = energy_cas_dyall(state_target) - energies(state_target)
|
||||
one_anhil_one_creat(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
@ -236,7 +229,7 @@ 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)]
|
||||
BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_act_orb,2,2,2,N_states)]
|
||||
implicit none
|
||||
integer :: i,j
|
||||
integer :: ispin,jspin,kspin
|
||||
@ -251,7 +244,6 @@ BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_a
|
||||
integer :: iorb,jorb
|
||||
integer :: korb
|
||||
integer :: state_target
|
||||
state_target = 1
|
||||
double precision :: energies(n_states_diag)
|
||||
do iorb = 1,n_act_orb
|
||||
do ispin = 1,2
|
||||
@ -278,6 +270,7 @@ BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_a
|
||||
enddo
|
||||
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)
|
||||
call apply_exc_to_psi(orb_k,hole_particle_k,spin_exc_k, &
|
||||
@ -285,7 +278,8 @@ BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_a
|
||||
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)
|
||||
two_anhil_one_creat(iorb,jorb,korb,ispin,jspin,kspin) = energy_cas_dyall(state_target) - energies(state_target)
|
||||
two_anhil_one_creat(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
@ -295,7 +289,7 @@ 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)]
|
||||
BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_act_orb,2,2,2,N_states)]
|
||||
implicit none
|
||||
integer :: i,j
|
||||
integer :: ispin,jspin,kspin
|
||||
@ -310,7 +304,6 @@ BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_a
|
||||
integer :: iorb,jorb
|
||||
integer :: korb
|
||||
integer :: state_target
|
||||
state_target = 1
|
||||
double precision :: energies(n_states_diag)
|
||||
do iorb = 1,n_act_orb
|
||||
do ispin = 1,2
|
||||
@ -336,6 +329,7 @@ BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_a
|
||||
psi_in_out(j,2,i) = psi_det(j,2,i)
|
||||
enddo
|
||||
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)
|
||||
call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, &
|
||||
@ -343,7 +337,8 @@ BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_a
|
||||
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)
|
||||
two_creat_one_anhil(iorb,jorb,korb,ispin,jspin,kspin) = energy_cas_dyall(state_target) - energies(state_target)
|
||||
two_creat_one_anhil(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
@ -353,7 +348,7 @@ BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_a
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2,2,2)]
|
||||
BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2,2,2,N_states)]
|
||||
implicit none
|
||||
integer :: i,j
|
||||
integer :: ispin,jspin,kspin
|
||||
@ -368,7 +363,6 @@ BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2
|
||||
integer :: iorb,jorb
|
||||
integer :: korb
|
||||
integer :: state_target
|
||||
state_target = 1
|
||||
double precision :: energies(n_states_diag)
|
||||
do iorb = 1,n_act_orb
|
||||
do ispin = 1,2
|
||||
@ -394,6 +388,7 @@ BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2
|
||||
psi_in_out(j,2,i) = psi_det(j,2,i)
|
||||
enddo
|
||||
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)
|
||||
call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, &
|
||||
@ -401,7 +396,8 @@ BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2
|
||||
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)
|
||||
three_creat(iorb,jorb,korb,ispin,jspin,kspin) = energy_cas_dyall(state_target) - energies(state_target)
|
||||
three_creat(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
@ -411,7 +407,7 @@ BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, three_anhil, (n_act_orb,n_act_orb,n_act_orb,2,2,2)]
|
||||
BEGIN_PROVIDER [ double precision, three_anhil, (n_act_orb,n_act_orb,n_act_orb,2,2,2,N_states)]
|
||||
implicit none
|
||||
integer :: i,j
|
||||
integer :: ispin,jspin,kspin
|
||||
@ -426,7 +422,6 @@ BEGIN_PROVIDER [ double precision, three_anhil, (n_act_orb,n_act_orb,n_act_orb,2
|
||||
integer :: iorb,jorb
|
||||
integer :: korb
|
||||
integer :: state_target
|
||||
state_target = 1
|
||||
double precision :: energies(n_states_diag)
|
||||
do iorb = 1,n_act_orb
|
||||
do ispin = 1,2
|
||||
@ -452,6 +447,7 @@ BEGIN_PROVIDER [ double precision, three_anhil, (n_act_orb,n_act_orb,n_act_orb,2
|
||||
psi_in_out(j,2,i) = psi_det(j,2,i)
|
||||
enddo
|
||||
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)
|
||||
call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, &
|
||||
@ -459,7 +455,8 @@ BEGIN_PROVIDER [ double precision, three_anhil, (n_act_orb,n_act_orb,n_act_orb,2
|
||||
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)
|
||||
three_anhil(iorb,jorb,korb,ispin,jspin,kspin) = energy_cas_dyall(state_target) - energies(state_target)
|
||||
three_anhil(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
@ -1,190 +0,0 @@
|
||||
|
||||
BEGIN_PROVIDER [ double precision, one_creat_spin_trace, (n_act_orb)]
|
||||
implicit none
|
||||
integer :: i
|
||||
do i = 1, n_act_orb
|
||||
one_creat_spin_trace(i) = one_creat(i,1) + one_creat(i,2)
|
||||
one_creat_spin_trace(i) = 0.5d0 * one_creat_spin_trace(i)
|
||||
enddo
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ double precision, one_anhil_spin_trace, (n_act_orb)]
|
||||
implicit none
|
||||
integer :: i
|
||||
do i = 1, n_act_orb
|
||||
one_anhil_spin_trace(i) = one_anhil(i,1) + one_anhil(i,2)
|
||||
one_anhil_spin_trace(i) = 0.5d0 * one_anhil_spin_trace(i)
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, two_creat_spin_trace, (n_act_orb,n_act_orb)]
|
||||
implicit none
|
||||
integer :: i,j
|
||||
integer :: ispin,jspin
|
||||
double precision :: counting
|
||||
do i = 1, n_act_orb
|
||||
do j = 1, n_act_orb
|
||||
two_creat_spin_trace(j,i) = 0.d0
|
||||
counting = 0.d0
|
||||
do ispin = 1, 2
|
||||
do jspin = 1,2
|
||||
two_creat_spin_trace(j,i) += two_creat(j,i,ispin,jspin)
|
||||
counting += 1.d0
|
||||
enddo
|
||||
enddo
|
||||
two_creat_spin_trace(j,i) = two_creat_spin_trace(j,i) / counting
|
||||
enddo
|
||||
enddo
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, two_anhil_spin_trace, (n_act_orb,n_act_orb)]
|
||||
implicit none
|
||||
integer :: i,j
|
||||
integer :: ispin,jspin
|
||||
double precision :: counting
|
||||
do i = 1, n_act_orb
|
||||
do j = 1, n_act_orb
|
||||
two_anhil_spin_trace(j,i) = 0.d0
|
||||
counting = 0.d0
|
||||
do ispin = 1, 2
|
||||
do jspin = 1,2
|
||||
two_anhil_spin_trace(j,i) += two_anhil(j,i,ispin,jspin)
|
||||
counting += 1.d0
|
||||
enddo
|
||||
enddo
|
||||
two_anhil_spin_trace(j,i) = two_anhil_spin_trace(j,i) / counting
|
||||
enddo
|
||||
enddo
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ double precision, one_anhil_one_creat_spin_trace, (n_act_orb,n_act_orb)]
|
||||
implicit none
|
||||
integer :: i,j
|
||||
integer :: ispin,jspin
|
||||
double precision :: counting
|
||||
do i = 1, n_act_orb
|
||||
do j = 1, n_act_orb
|
||||
one_anhil_one_creat_spin_trace(j,i) = 0.d0
|
||||
counting = 0.d0
|
||||
do ispin = 1, 2
|
||||
do jspin = 1,2
|
||||
one_anhil_one_creat_spin_trace(j,i) += one_anhil_one_creat(j,i,jspin,ispin)
|
||||
counting += 1.d0
|
||||
enddo
|
||||
enddo
|
||||
one_anhil_one_creat_spin_trace(j,i) = one_anhil_one_creat_spin_trace(j,i) / counting
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ double precision, two_anhil_one_creat_spin_trace, (n_act_orb,n_act_orb,n_act_orb)]
|
||||
implicit none
|
||||
integer :: i,j,k
|
||||
integer :: ispin,jspin,kspin
|
||||
double precision :: counting
|
||||
|
||||
do i = 1, n_act_orb
|
||||
do j = 1, n_act_orb
|
||||
do k = 1, n_act_orb
|
||||
two_anhil_one_creat_spin_trace(k,j,i) = 0.d0
|
||||
counting = 0.d0
|
||||
do ispin = 1, 2
|
||||
do jspin = 1,2
|
||||
do kspin = 1,2
|
||||
two_anhil_one_creat_spin_trace(k,j,i) += two_anhil_one_creat(k,j,i,kspin,jspin,ispin)
|
||||
counting += 1.d0
|
||||
enddo
|
||||
enddo
|
||||
two_anhil_one_creat_spin_trace(k,j,i) = two_anhil_one_creat_spin_trace(k,j,i) / counting
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, two_creat_one_anhil_spin_trace, (n_act_orb,n_act_orb,n_act_orb)]
|
||||
implicit none
|
||||
integer :: i,j,k
|
||||
integer :: ispin,jspin,kspin
|
||||
double precision :: counting
|
||||
|
||||
do i = 1, n_act_orb
|
||||
do j = 1, n_act_orb
|
||||
do k = 1, n_act_orb
|
||||
two_creat_one_anhil_spin_trace(k,j,i) = 0.d0
|
||||
counting = 0.d0
|
||||
do ispin = 1, 2
|
||||
do jspin = 1,2
|
||||
do kspin = 1,2
|
||||
two_creat_one_anhil_spin_trace(k,j,i) += two_creat_one_anhil(k,j,i,kspin,jspin,ispin)
|
||||
counting += 1.d0
|
||||
enddo
|
||||
enddo
|
||||
two_creat_one_anhil_spin_trace(k,j,i) = two_creat_one_anhil_spin_trace(k,j,i) / counting
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ double precision, three_creat_spin_trace, (n_act_orb,n_act_orb,n_act_orb)]
|
||||
implicit none
|
||||
integer :: i,j,k
|
||||
integer :: ispin,jspin,kspin
|
||||
double precision :: counting
|
||||
|
||||
do i = 1, n_act_orb
|
||||
do j = 1, n_act_orb
|
||||
do k = 1, n_act_orb
|
||||
three_creat_spin_trace(k,j,i) = 0.d0
|
||||
counting = 0.d0
|
||||
do ispin = 1, 2
|
||||
do jspin = 1,2
|
||||
do kspin = 1,2
|
||||
three_creat_spin_trace(k,j,i) += three_creat(k,j,i,kspin,jspin,ispin)
|
||||
counting += 1.d0
|
||||
enddo
|
||||
enddo
|
||||
three_creat_spin_trace(k,j,i) = three_creat_spin_trace(k,j,i) / counting
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ double precision, three_anhil_spin_trace, (n_act_orb,n_act_orb,n_act_orb)]
|
||||
implicit none
|
||||
integer :: i,j,k
|
||||
integer :: ispin,jspin,kspin
|
||||
double precision :: counting
|
||||
|
||||
do i = 1, n_act_orb
|
||||
do j = 1, n_act_orb
|
||||
do k = 1, n_act_orb
|
||||
three_anhil_spin_trace(k,j,i) = 0.d0
|
||||
counting = 0.d0
|
||||
do ispin = 1, 2
|
||||
do jspin = 1,2
|
||||
do kspin = 1,2
|
||||
three_anhil_spin_trace(k,j,i) += three_anhil(k,j,i,kspin,jspin,ispin)
|
||||
counting += 1.d0
|
||||
enddo
|
||||
enddo
|
||||
three_anhil_spin_trace(k,j,i) = three_anhil_spin_trace(k,j,i) / counting
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
@ -44,7 +44,7 @@
|
||||
enddo
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [double precision, fock_core_inactive_from_act, (mo_tot_num,2)]
|
||||
BEGIN_PROVIDER [double precision, fock_core_inactive_from_act, (mo_tot_num,2,N_states)]
|
||||
BEGIN_DOC
|
||||
! inactive part of the fock operator with contributions only from the active
|
||||
END_DOC
|
||||
@ -55,15 +55,17 @@
|
||||
double precision :: coulomb, exchange
|
||||
double precision :: get_mo_bielec_integral_schwartz
|
||||
integer :: j_act_orb,k_act_orb,i_inact_core_orb
|
||||
integer :: i_state
|
||||
|
||||
do i_state = 1,N_states
|
||||
do i = 1, n_core_inact_orb
|
||||
accu_coulomb = 0.d0
|
||||
accu_exchange = 0.d0
|
||||
i_inact_core_orb = list_core_inact(i)
|
||||
do j = 1, n_act_orb
|
||||
j_act_orb = list_act(j)
|
||||
na = one_body_dm_mo_alpha(j_act_orb,j_act_orb)
|
||||
nb = one_body_dm_mo_beta(j_act_orb,j_act_orb)
|
||||
na = one_body_dm_mo_alpha(j_act_orb,j_act_orb,i_state)
|
||||
nb = one_body_dm_mo_beta(j_act_orb,j_act_orb,i_state)
|
||||
ntot = na + nb
|
||||
coulomb = mo_bielec_integral_jj(i_inact_core_orb,j_act_orb)
|
||||
exchange = mo_bielec_integral_jj_exchange(i_inact_core_orb,j_act_orb)
|
||||
@ -72,8 +74,8 @@
|
||||
accu_exchange(2) += nb * exchange
|
||||
do k = j+1, n_act_orb
|
||||
k_act_orb = list_act(k)
|
||||
na = one_body_dm_mo_alpha(j_act_orb,k_act_orb)
|
||||
nb = one_body_dm_mo_beta(j_act_orb,k_act_orb)
|
||||
na = one_body_dm_mo_alpha(j_act_orb,k_act_orb,i_state)
|
||||
nb = one_body_dm_mo_beta(j_act_orb,k_act_orb,i_state)
|
||||
ntot = na + nb
|
||||
coulomb = get_mo_bielec_integral_schwartz(j_act_orb,i_inact_core_orb,k_act_orb,i_inact_core_orb,mo_integrals_map)
|
||||
exchange = get_mo_bielec_integral_schwartz(j_act_orb,k_act_orb,i_inact_core_orb,i_inact_core_orb,mo_integrals_map)
|
||||
@ -82,12 +84,13 @@
|
||||
accu_exchange(2) += 2.d0 * nb * exchange
|
||||
enddo
|
||||
enddo
|
||||
fock_core_inactive_from_act(i_inact_core_orb,1) = accu_coulomb + accu_exchange(1)
|
||||
fock_core_inactive_from_act(i_inact_core_orb,2) = accu_coulomb + accu_exchange(2)
|
||||
fock_core_inactive_from_act(i_inact_core_orb,1,i_state) = accu_coulomb + accu_exchange(1)
|
||||
fock_core_inactive_from_act(i_inact_core_orb,2,i_state) = accu_coulomb + accu_exchange(2)
|
||||
enddo
|
||||
enddo
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [double precision, fock_virt_from_act, (mo_tot_num,2)]
|
||||
BEGIN_PROVIDER [double precision, fock_virt_from_act, (mo_tot_num,2,N_states)]
|
||||
BEGIN_DOC
|
||||
! virtual part of the fock operator with contributions only from the active
|
||||
END_DOC
|
||||
@ -98,15 +101,18 @@
|
||||
double precision :: coulomb, exchange
|
||||
double precision :: get_mo_bielec_integral_schwartz
|
||||
integer :: j_act_orb,i_virt_orb,k_act_orb
|
||||
integer :: i_state
|
||||
! TODO : inverse loop of i_state
|
||||
|
||||
do i_state = 1, N_states
|
||||
do i = 1, n_virt_orb
|
||||
accu_coulomb = 0.d0
|
||||
accu_exchange = 0.d0
|
||||
i_virt_orb = list_virt(i)
|
||||
do j = 1, n_act_orb
|
||||
j_act_orb = list_act(j)
|
||||
na = one_body_dm_mo_alpha(j_act_orb,j_act_orb)
|
||||
nb = one_body_dm_mo_beta(j_act_orb,j_act_orb)
|
||||
na = one_body_dm_mo_alpha(j_act_orb,j_act_orb,i_state)
|
||||
nb = one_body_dm_mo_beta(j_act_orb,j_act_orb,i_state)
|
||||
ntot = na + nb
|
||||
coulomb = mo_bielec_integral_jj(i_virt_orb,j_act_orb)
|
||||
exchange = mo_bielec_integral_jj_exchange(i_virt_orb,j_act_orb)
|
||||
@ -115,8 +121,8 @@
|
||||
accu_exchange(2) += nb * exchange
|
||||
do k = j+1, n_act_orb
|
||||
k_act_orb = list_act(k)
|
||||
na = one_body_dm_mo_alpha(j_act_orb,k_act_orb)
|
||||
nb = one_body_dm_mo_beta(j_act_orb,k_act_orb)
|
||||
na = one_body_dm_mo_alpha(j_act_orb,k_act_orb,i_state)
|
||||
nb = one_body_dm_mo_beta(j_act_orb,k_act_orb,i_state)
|
||||
ntot = na + nb
|
||||
coulomb = get_mo_bielec_integral_schwartz(j_act_orb,i_virt_orb,k_act_orb,i_virt_orb,mo_integrals_map)
|
||||
exchange = get_mo_bielec_integral_schwartz(j_act_orb,k_act_orb,i_virt_orb,i_virt_orb,mo_integrals_map)
|
||||
@ -125,40 +131,47 @@
|
||||
accu_exchange(2) += 2.d0 * nb * exchange
|
||||
enddo
|
||||
enddo
|
||||
fock_virt_from_act(i_virt_orb,1) = accu_coulomb + accu_exchange(1)
|
||||
fock_virt_from_act(i_virt_orb,2) = accu_coulomb + accu_exchange(2)
|
||||
fock_virt_from_act(i_virt_orb,1,i_state) = accu_coulomb + accu_exchange(1)
|
||||
fock_virt_from_act(i_virt_orb,2,i_state) = accu_coulomb + accu_exchange(2)
|
||||
enddo
|
||||
enddo
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [double precision, fock_core_inactive_total, (mo_tot_num,2)]
|
||||
&BEGIN_PROVIDER [double precision, fock_core_inactive_total_spin_trace, (mo_tot_num)]
|
||||
BEGIN_PROVIDER [double precision, fock_core_inactive_total, (mo_tot_num,2,N_states)]
|
||||
&BEGIN_PROVIDER [double precision, fock_core_inactive_total_spin_trace, (mo_tot_num,N_states)]
|
||||
BEGIN_DOC
|
||||
! inactive part of the fock operator
|
||||
END_DOC
|
||||
implicit none
|
||||
integer :: i
|
||||
integer :: i_inact_core_orb
|
||||
integer :: i_state
|
||||
do i_state = 1, N_states
|
||||
do i = 1, n_core_inact_orb
|
||||
i_inact_core_orb = list_core_inact(i)
|
||||
fock_core_inactive_total(i_inact_core_orb,1) = fock_core_inactive(i_inact_core_orb) + fock_core_inactive_from_act(i_inact_core_orb,1)
|
||||
fock_core_inactive_total(i_inact_core_orb,2) = fock_core_inactive(i_inact_core_orb) + fock_core_inactive_from_act(i_inact_core_orb,2)
|
||||
fock_core_inactive_total_spin_trace(i_inact_core_orb) = 0.5d0 * (fock_core_inactive_total(i_inact_core_orb,1) + fock_core_inactive_total(i_inact_core_orb,2))
|
||||
fock_core_inactive_total(i_inact_core_orb,1,i_state) = fock_core_inactive(i_inact_core_orb) + fock_core_inactive_from_act(i_inact_core_orb,1,i_state)
|
||||
fock_core_inactive_total(i_inact_core_orb,2,i_state) = fock_core_inactive(i_inact_core_orb) + fock_core_inactive_from_act(i_inact_core_orb,2,i_state)
|
||||
fock_core_inactive_total_spin_trace(i_inact_core_orb,i_state) = 0.5d0 * (fock_core_inactive_total(i_inact_core_orb,1,i_state) + fock_core_inactive_total(i_inact_core_orb,2,i_state))
|
||||
enddo
|
||||
enddo
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [double precision, fock_virt_total, (mo_tot_num,2)]
|
||||
&BEGIN_PROVIDER [double precision, fock_virt_total_spin_trace, (mo_tot_num)]
|
||||
BEGIN_PROVIDER [double precision, fock_virt_total, (mo_tot_num,2,N_states)]
|
||||
&BEGIN_PROVIDER [double precision, fock_virt_total_spin_trace, (mo_tot_num,N_states)]
|
||||
BEGIN_DOC
|
||||
! inactive part of the fock operator
|
||||
END_DOC
|
||||
implicit none
|
||||
integer :: i
|
||||
integer :: i_virt_orb
|
||||
integer :: i_state
|
||||
do i_state = 1, N_states
|
||||
do i = 1, n_virt_orb
|
||||
i_virt_orb= list_virt(i)
|
||||
fock_virt_total(i_virt_orb,1) = fock_virt_from_core_inact(i_virt_orb) + fock_virt_from_act(i_virt_orb,1)+ mo_mono_elec_integral(i_virt_orb,i_virt_orb)
|
||||
fock_virt_total(i_virt_orb,2) = fock_virt_from_core_inact(i_virt_orb) + fock_virt_from_act(i_virt_orb,2)+ mo_mono_elec_integral(i_virt_orb,i_virt_orb)
|
||||
fock_virt_total_spin_trace(i_virt_orb) = 0.5d0 * ( fock_virt_total(i_virt_orb,1) + fock_virt_total(i_virt_orb,2) )
|
||||
fock_virt_total(i_virt_orb,1,i_state) = fock_virt_from_core_inact(i_virt_orb) + fock_virt_from_act(i_virt_orb,1,i_state)+ mo_mono_elec_integral(i_virt_orb,i_virt_orb)
|
||||
fock_virt_total(i_virt_orb,2,i_state) = fock_virt_from_core_inact(i_virt_orb) + fock_virt_from_act(i_virt_orb,2,i_state)+ mo_mono_elec_integral(i_virt_orb,i_virt_orb)
|
||||
fock_virt_total_spin_trace(i_virt_orb,i_state) = 0.5d0 * ( fock_virt_total(i_virt_orb,1,i_state) + fock_virt_total(i_virt_orb,2,i_state) )
|
||||
enddo
|
||||
enddo
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
@ -30,19 +30,19 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip
|
||||
integer :: degree_alpha(psi_det_size)
|
||||
logical :: fullMatch
|
||||
|
||||
double precision :: delta_e_array(psi_det_size)
|
||||
double precision :: delta_e_inv_array(psi_det_size,N_states)
|
||||
double precision :: hij_array(psi_det_size)
|
||||
|
||||
integer(bit_kind) :: tq(Nint,2,n_selected)
|
||||
integer :: N_tq
|
||||
|
||||
double precision :: hialpha
|
||||
double precision :: hialpha,hij
|
||||
integer :: i_state, i_alpha
|
||||
|
||||
integer(bit_kind),allocatable :: miniList(:,:,:)
|
||||
integer,allocatable :: idx_miniList(:)
|
||||
integer :: N_miniList, leng
|
||||
double precision :: delta_e_final,hij_tmp
|
||||
double precision :: delta_e(N_states),hij_tmp
|
||||
integer :: index_i,index_j
|
||||
|
||||
|
||||
@ -76,12 +76,12 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip
|
||||
! coef_pert = 0.d0
|
||||
do i = 1,idx_alpha(0)
|
||||
index_i = idx_alpha(i)
|
||||
call get_delta_e_dyall(psi_det(1,1,index_i),tq(1,1,i_alpha),delta_e_final)
|
||||
call get_delta_e_dyall(psi_det(1,1,index_i),tq(1,1,i_alpha),delta_e)
|
||||
call i_h_j(tq(1,1,i_alpha),psi_det(1,1,index_i),Nint,hialpha)
|
||||
delta_e_array(index_i) = 1.d0/delta_e_final
|
||||
hij_array(index_i) = hialpha
|
||||
! ihpsi0 += hialpha * psi_coef(index_i,1)
|
||||
! coef_pert += hialpha * psi_coef(index_i,1) * delta_e_array(index_i)
|
||||
do i_state = 1,N_states
|
||||
delta_e_inv_array(index_i,i_state) = 1.d0/delta_e(i_state)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
do i=1,idx_alpha(0)
|
||||
@ -91,7 +91,8 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip
|
||||
do j = 1, idx_alpha(0)
|
||||
index_j = idx_alpha(j)
|
||||
do i_state=1,N_states
|
||||
delta_ij_(index_i,index_j,i_state) += hij_array(index_j) * hij_tmp * delta_e_array(index_j)
|
||||
! standard dressing first order
|
||||
delta_ij_(index_i,index_j,i_state) += hij_array(index_j) * hij_tmp * delta_e_inv_array(index_j,i_state)
|
||||
enddo
|
||||
enddo
|
||||
call omp_unset_lock( psi_ref_bis_lock(index_i))
|
||||
|
||||
@ -14,7 +14,8 @@
|
||||
! Dressing matrix in N_det basis
|
||||
END_DOC
|
||||
integer :: i,j,m
|
||||
double precision :: accu
|
||||
integer :: i_state
|
||||
double precision :: accu(N_states)
|
||||
double precision, allocatable :: delta_ij_tmp(:,:,:)
|
||||
|
||||
|
||||
@ -27,118 +28,136 @@
|
||||
delta_ij_tmp = 0.d0
|
||||
call H_apply_mrpt_1h(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 += delta_ij_tmp(j,i,1) * psi_coef(i,1) * psi_coef(j,1)
|
||||
delta_ij(j,i,1) += delta_ij_tmp(j,i,1)
|
||||
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
|
||||
enddo
|
||||
second_order_pt_new_1h(i_state) = accu(i_state)
|
||||
enddo
|
||||
print*, '1h = ',accu
|
||||
second_order_pt_new_1h(1) = 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 += delta_ij_tmp(j,i,1) * psi_coef(i,1) * psi_coef(j,1)
|
||||
delta_ij(j,i,1) += delta_ij_tmp(j,i,1)
|
||||
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
|
||||
enddo
|
||||
second_order_pt_new_1p(i_state) = accu(i_state)
|
||||
enddo
|
||||
print*, '1p = ',accu
|
||||
second_order_pt_new_1p(1) = accu
|
||||
|
||||
! 1h1p
|
||||
delta_ij_tmp = 0.d0
|
||||
call H_apply_mrpt_1h1p(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 += delta_ij_tmp(j,i,1) * psi_coef(i,1) * psi_coef(j,1)
|
||||
delta_ij(j,i,1) += delta_ij_tmp(j,i,1)
|
||||
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
|
||||
enddo
|
||||
second_order_pt_new_1h1p(i_state) = accu(i_state)
|
||||
enddo
|
||||
print*, '1h1p = ',accu
|
||||
second_order_pt_new_1h1p(1) = 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 += delta_ij_tmp(j,i,1) * psi_coef(i,1) * psi_coef(j,1)
|
||||
delta_ij(j,i,1) += delta_ij_tmp(j,i,1)
|
||||
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
|
||||
enddo
|
||||
second_order_pt_new_2h(i_state) = accu(i_state)
|
||||
enddo
|
||||
print*, '2h = ',accu
|
||||
second_order_pt_new_2h(1) = 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 += delta_ij_tmp(j,i,1) * psi_coef(i,1) * psi_coef(j,1)
|
||||
delta_ij(j,i,1) += delta_ij_tmp(j,i,1)
|
||||
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
|
||||
enddo
|
||||
second_order_pt_new_2p(i_state) = accu(i_state)
|
||||
enddo
|
||||
print*, '2p = ',accu
|
||||
second_order_pt_new_2p(1) = accu
|
||||
|
||||
! 1h2p
|
||||
delta_ij_tmp = 0.d0
|
||||
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 += delta_ij_tmp(j,i,1) * psi_coef(i,1) * psi_coef(j,1)
|
||||
delta_ij(j,i,1) += delta_ij_tmp(j,i,1)
|
||||
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
|
||||
enddo
|
||||
second_order_pt_new_1h2p(i_state) = accu(i_state)
|
||||
enddo
|
||||
print*, '1h2p = ',accu
|
||||
second_order_pt_new_1h2p(1) = accu
|
||||
|
||||
! 2h1p
|
||||
delta_ij_tmp = 0.d0
|
||||
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 += delta_ij_tmp(j,i,1) * psi_coef(i,1) * psi_coef(j,1)
|
||||
delta_ij(j,i,1) += delta_ij_tmp(j,i,1)
|
||||
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
|
||||
enddo
|
||||
second_order_pt_new_2h1p(i_state) = accu(i_state)
|
||||
enddo
|
||||
print*, '2h1p = ',accu
|
||||
second_order_pt_new_2h1p(1) = accu
|
||||
|
||||
! 2h2p
|
||||
delta_ij_tmp = 0.d0
|
||||
call H_apply_mrpt_2h2p(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 += delta_ij_tmp(j,i,1) * psi_coef(i,1) * psi_coef(j,1)
|
||||
delta_ij(j,i,1) += delta_ij_tmp(j,i,1)
|
||||
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
|
||||
enddo
|
||||
second_order_pt_new_2h2p(i_state) = accu(i_state)
|
||||
enddo
|
||||
print*, '2h2p = ',accu
|
||||
second_order_pt_new_2h2p(1) = accu
|
||||
|
||||
! total
|
||||
accu = 0.d0
|
||||
do i_state = 1, N_states
|
||||
do i = 1, N_det
|
||||
do j = 1, N_det
|
||||
accu += delta_ij(j,i,1) * psi_coef(i,1) * psi_coef(j,1)
|
||||
write(*,'(1000(F16.10,x))')delta_ij(i,:,:)
|
||||
do j = i_state, N_det
|
||||
accu(i_state) += delta_ij(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
|
||||
enddo
|
||||
enddo
|
||||
print*, 'total= ',accu
|
||||
second_order_pt_new(1) = accu
|
||||
second_order_pt_new(i_state) = accu(i_state)
|
||||
print*, 'total= ',accu(i_state)
|
||||
enddo
|
||||
|
||||
|
||||
! write(*,'(1000(F16.10,x))')delta_ij(i,:,:)
|
||||
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
@ -152,33 +152,51 @@ subroutine give_particles_in_virt_space(det_1,n_particles_spin,n_particles,parti
|
||||
end
|
||||
|
||||
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
|
||||
!
|
||||
! 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 from the contracted application of the excitation
|
||||
!
|
||||
! operator in the active space that lead from det_1 to det_2
|
||||
END_DOC
|
||||
implicit none
|
||||
use bitmasks
|
||||
double precision, intent(out) :: delta_e_final
|
||||
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
|
||||
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)
|
||||
delta_e_inactive += fock_core_inactive_total_spin_trace(i_hole_inact)
|
||||
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)
|
||||
delta_e_inactive += fock_core_inactive_total_spin_trace(i_hole_inact)
|
||||
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
|
||||
double precision :: delta_e_virt(N_states)
|
||||
integer :: i_part_virt
|
||||
integer :: n_particles_spin(2)
|
||||
integer :: n_particles
|
||||
@ -188,12 +206,16 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final)
|
||||
delta_e_virt = 0.d0
|
||||
do i = 1, n_particles_spin(1)
|
||||
i_part_virt = particles_list(i,1)
|
||||
delta_e_virt += fock_virt_total_spin_trace(i_part_virt)
|
||||
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)
|
||||
delta_e_virt += fock_virt_total_spin_trace(i_part_virt)
|
||||
do i_state = 1, N_states
|
||||
delta_e_virt += fock_virt_total_spin_trace(i_part_virt,i_state)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
|
||||
@ -203,7 +225,7 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final)
|
||||
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
|
||||
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)
|
||||
@ -265,14 +287,18 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final)
|
||||
! delta_e_act += one_creat_spin_trace(i_particle_act )
|
||||
ispin = particle_list_practical(1,1)
|
||||
i_particle_act = particle_list_practical(2,1)
|
||||
delta_e_act += one_creat(i_particle_act,ispin)
|
||||
do i_state = 1, N_states
|
||||
delta_e_act(i_state) += one_creat(i_particle_act,ispin,i_state)
|
||||
enddo
|
||||
|
||||
else if (n_holes_act == 1 .and. n_particles_act == 0) then
|
||||
! i_hole_act = holes_active_list_spin_traced(1)
|
||||
! delta_e_act += one_anhil_spin_trace(i_hole_act )
|
||||
ispin = hole_list_practical(1,1)
|
||||
i_hole_act = hole_list_practical(2,1)
|
||||
delta_e_act += one_anhil(i_hole_act , ispin)
|
||||
do i_state = 1, N_states
|
||||
delta_e_act(i_state) += one_anhil(i_hole_act , ispin,i_state)
|
||||
enddo
|
||||
|
||||
else if (n_holes_act == 1 .and. n_particles_act == 1) then
|
||||
! i_hole_act = holes_active_list_spin_traced(1)
|
||||
@ -284,7 +310,9 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final)
|
||||
! first particle
|
||||
jspin = particle_list_practical(1,1)
|
||||
i_particle_act = particle_list_practical(2,1)
|
||||
delta_e_act += one_anhil_one_creat(i_particle_act,i_hole_act,jspin,ispin)
|
||||
do i_state = 1, N_states
|
||||
delta_e_act(i_state) += one_anhil_one_creat(i_particle_act,i_hole_act,jspin,ispin,i_state)
|
||||
enddo
|
||||
|
||||
else if (n_holes_act == 2 .and. n_particles_act == 0) then
|
||||
! i_hole_act = holes_active_list_spin_traced(1)
|
||||
@ -294,7 +322,9 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final)
|
||||
i_hole_act = hole_list_practical(2,1)
|
||||
jspin = hole_list_practical(1,2)
|
||||
j_hole_act = hole_list_practical(2,2)
|
||||
delta_e_act += two_anhil(i_hole_act,j_hole_act,ispin,jspin)
|
||||
do i_state = 1, N_states
|
||||
delta_e_act(i_state) += two_anhil(i_hole_act,j_hole_act,ispin,jspin,i_state)
|
||||
enddo
|
||||
|
||||
else if (n_holes_act == 0 .and. n_particles_act == 2) then
|
||||
! i_particle_act = particles_active_list_spin_traced(1)
|
||||
@ -304,7 +334,9 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final)
|
||||
i_particle_act = particle_list_practical(2,1)
|
||||
jspin = particle_list_practical(1,2)
|
||||
j_particle_act = particle_list_practical(2,2)
|
||||
delta_e_act += two_creat(i_particle_act,j_particle_act,ispin,jspin)
|
||||
do i_state = 1, N_states
|
||||
delta_e_act(i_state) += two_creat(i_particle_act,j_particle_act,ispin,jspin,i_state)
|
||||
enddo
|
||||
|
||||
else if (n_holes_act == 2 .and. n_particles_act == 1) then
|
||||
! i_hole_act = holes_active_list_spin_traced(1)
|
||||
@ -324,7 +356,9 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final)
|
||||
! first particle
|
||||
kspin = particle_list_practical(1,1)
|
||||
i_particle_act = particle_list_practical(2,1)
|
||||
delta_e_act += two_anhil_one_creat(i_particle_act,i_hole_act,j_hole_act,kspin,ispin,jspin)
|
||||
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)
|
||||
enddo
|
||||
|
||||
else if (n_holes_act == 1 .and. n_particles_act == 2) then
|
||||
! i_hole_act = holes_active_list_spin_traced(1)
|
||||
@ -342,7 +376,9 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final)
|
||||
kspin = particle_list_practical(1,2)
|
||||
j_particle_act = particle_list_practical(2,2)
|
||||
|
||||
delta_e_act += two_creat_one_anhil(i_particle_act,j_particle_act,i_hole_act,jspin,kspin,ispin)
|
||||
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)
|
||||
enddo
|
||||
|
||||
else if (n_holes_act == 3 .and. n_particles_act == 0) then
|
||||
! i_hole_act = holes_active_list_spin_traced(1)
|
||||
@ -359,7 +395,9 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final)
|
||||
! third hole
|
||||
kspin = hole_list_practical(1,3)
|
||||
k_hole_act = hole_list_practical(2,3)
|
||||
delta_e_act += three_anhil(i_hole_act,j_hole_act,k_hole_act,ispin,jspin,kspin)
|
||||
do i_state = 1, N_states
|
||||
delta_e_act(i_state) += three_anhil(i_hole_act,j_hole_act,k_hole_act,ispin,jspin,kspin,i_state)
|
||||
enddo
|
||||
|
||||
else if (n_holes_act == 0 .and. n_particles_act == 3) then
|
||||
! i_particle_act = particles_active_list_spin_traced(1)
|
||||
@ -376,7 +414,9 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final)
|
||||
kspin = particle_list_practical(1,3)
|
||||
k_particle_act = particle_list_practical(2,3)
|
||||
|
||||
delta_e_act += three_creat(i_particle_act,j_particle_act,k_particle_act,ispin,jspin,kspin)
|
||||
do i_state = 1, N_states
|
||||
delta_e_act(i_state) += three_creat(i_particle_act,j_particle_act,k_particle_act,ispin,jspin,kspin,i_state)
|
||||
enddo
|
||||
|
||||
else if (n_holes_act .ge. 2 .and. n_particles_act .ge.2) then
|
||||
|
||||
@ -388,265 +428,8 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final)
|
||||
!print*, one_anhil_spin_trace(1), one_anhil_spin_trace(2)
|
||||
|
||||
|
||||
delta_e_final = delta_e_act + delta_e_inactive - delta_e_virt
|
||||
!if(delta_e_final .le. -100d0.or.delta_e_final > 0.d0 .or. delta_e_final == 0.d0)then
|
||||
!if(delta_e_final == 0.d0)then
|
||||
if(.False.)then
|
||||
call debug_det(det_1,N_int)
|
||||
call debug_det(det_2,N_int)
|
||||
print*, 'n_holes_act,n_particles_act'
|
||||
print*, n_holes_act,n_particles_act
|
||||
print*, 'delta_e_act,delta_e_inactive,delta_e_vir'
|
||||
print*, delta_e_act,delta_e_inactive,delta_e_virt
|
||||
delta_e_final = -1000.d0
|
||||
!stop
|
||||
|
||||
endif
|
||||
|
||||
end
|
||||
|
||||
subroutine get_delta_e_dyall_verbose(det_1,det_2,delta_e_final)
|
||||
implicit none
|
||||
use bitmasks
|
||||
double precision, intent(out) :: delta_e_final
|
||||
integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2)
|
||||
integer :: i,j,k,l
|
||||
|
||||
integer :: n_holes_spin(2)
|
||||
integer :: n_holes
|
||||
integer :: holes_list(N_int*bit_kind_size,2)
|
||||
|
||||
|
||||
double precision :: delta_e_inactive
|
||||
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)
|
||||
delta_e_inactive += fock_core_inactive_total_spin_trace(i_hole_inact)
|
||||
enddo
|
||||
|
||||
do i = 1, n_holes_spin(2)
|
||||
i_hole_inact = holes_list(i,2)
|
||||
delta_e_inactive += fock_core_inactive_total_spin_trace(i_hole_inact)
|
||||
enddo
|
||||
|
||||
double precision :: delta_e_virt
|
||||
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)
|
||||
delta_e_virt += fock_virt_total_spin_trace(i_part_virt)
|
||||
enddo
|
||||
|
||||
do i = 1, n_particles_spin(2)
|
||||
i_part_virt = particles_list(i,2)
|
||||
delta_e_virt += fock_virt_total_spin_trace(i_part_virt)
|
||||
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
|
||||
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
|
||||
hole_list_practical(2,icountbis) = holes_active_list(i,1)
|
||||
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
|
||||
if (n_holes_act == 0 .and. n_particles_act == 1) then
|
||||
! i_particle_act = particles_active_list_spin_traced(1)
|
||||
! delta_e_act += one_creat_spin_trace(i_particle_act )
|
||||
ispin = particle_list_practical(1,1)
|
||||
i_particle_act = particle_list_practical(2,1)
|
||||
delta_e_act += one_creat(i_particle_act,ispin)
|
||||
|
||||
else if (n_holes_act == 1 .and. n_particles_act == 0) then
|
||||
! i_hole_act = holes_active_list_spin_traced(1)
|
||||
! delta_e_act += one_anhil_spin_trace(i_hole_act )
|
||||
ispin = hole_list_practical(1,1)
|
||||
i_hole_act = hole_list_practical(2,1)
|
||||
delta_e_act += one_anhil(i_hole_act , ispin)
|
||||
|
||||
else if (n_holes_act == 1 .and. n_particles_act == 1) then
|
||||
! i_hole_act = holes_active_list_spin_traced(1)
|
||||
! i_particle_act = particles_active_list_spin_traced(1)
|
||||
! delta_e_act += one_anhil_one_creat_spin_trace(i_hole_act,i_particle_act)
|
||||
! first hole
|
||||
ispin = hole_list_practical(1,1)
|
||||
i_hole_act = hole_list_practical(2,1)
|
||||
! first particle
|
||||
jspin = particle_list_practical(1,1)
|
||||
i_particle_act = particle_list_practical(2,1)
|
||||
delta_e_act += one_anhil_one_creat(i_particle_act,i_hole_act,jspin,ispin)
|
||||
|
||||
else if (n_holes_act == 2 .and. n_particles_act == 0) then
|
||||
! i_hole_act = holes_active_list_spin_traced(1)
|
||||
! j_hole_act = holes_active_list_spin_traced(1)
|
||||
! delta_e_act += two_anhil_spin_trace(i_hole_act,j_hole_act)
|
||||
ispin = hole_list_practical(1,1)
|
||||
i_hole_act = hole_list_practical(2,1)
|
||||
jspin = hole_list_practical(1,2)
|
||||
j_hole_act = hole_list_practical(2,2)
|
||||
delta_e_act += two_anhil(i_hole_act,j_hole_act,ispin,jspin)
|
||||
|
||||
else if (n_holes_act == 0 .and. n_particles_act == 2) then
|
||||
! i_particle_act = particles_active_list_spin_traced(1)
|
||||
! j_particle_act = particles_active_list_spin_traced(2)
|
||||
! delta_e_act += two_creat_spin_trace(i_particle_act,j_particle_act)
|
||||
ispin = particle_list_practical(1,1)
|
||||
i_particle_act = particle_list_practical(2,1)
|
||||
jspin = particle_list_practical(1,2)
|
||||
j_particle_act = particle_list_practical(2,2)
|
||||
delta_e_act += two_creat(i_particle_act,j_particle_act,ispin,jspin)
|
||||
|
||||
else if (n_holes_act == 2 .and. n_particles_act == 1) then
|
||||
! i_hole_act = holes_active_list_spin_traced(1)
|
||||
! j_hole_act = holes_active_list_spin_traced(2)
|
||||
! i_particle_act = particles_active_list_spin_traced(1)
|
||||
! print*, 'i_hole_act,j_hole_act,i_particle_act'
|
||||
! print*, i_hole_act,j_hole_act,i_particle_act
|
||||
! print*, two_anhil_one_creat_spin_trace(i_hole_act,j_hole_act,i_particle_act)
|
||||
! delta_e_act += two_anhil_one_creat_spin_trace(i_hole_act,j_hole_act,i_particle_act)
|
||||
|
||||
! first hole
|
||||
ispin = hole_list_practical(1,1)
|
||||
i_hole_act = hole_list_practical(2,1)
|
||||
! second hole
|
||||
jspin = hole_list_practical(1,2)
|
||||
j_hole_act = hole_list_practical(2,2)
|
||||
! first particle
|
||||
kspin = particle_list_practical(1,1)
|
||||
i_particle_act = particle_list_practical(2,1)
|
||||
delta_e_act += two_anhil_one_creat(i_particle_act,i_hole_act,j_hole_act,kspin,ispin,jspin)
|
||||
|
||||
else if (n_holes_act == 1 .and. n_particles_act == 2) then
|
||||
! i_hole_act = holes_active_list_spin_traced(1)
|
||||
! i_particle_act = particles_active_list_spin_traced(1)
|
||||
! j_particle_act = particles_active_list_spin_traced(2)
|
||||
! delta_e_act += two_creat_one_anhil_spin_trace(i_hole_act,i_particle_act,j_particle_act)
|
||||
|
||||
! first hole
|
||||
ispin = hole_list_practical(1,1)
|
||||
i_hole_act = hole_list_practical(2,1)
|
||||
! first particle
|
||||
jspin = particle_list_practical(1,1)
|
||||
i_particle_act = particle_list_practical(2,1)
|
||||
! second particle
|
||||
kspin = particle_list_practical(1,2)
|
||||
j_particle_act = particle_list_practical(2,2)
|
||||
|
||||
delta_e_act += two_creat_one_anhil(i_particle_act,j_particle_act,i_hole_act,jspin,kspin,ispin)
|
||||
|
||||
else if (n_holes_act == 3 .and. n_particles_act == 0) then
|
||||
! i_hole_act = holes_active_list_spin_traced(1)
|
||||
! j_hole_act = holes_active_list_spin_traced(2)
|
||||
! k_hole_act = holes_active_list_spin_traced(3)
|
||||
! delta_e_act += three_anhil_spin_trace(i_hole_act,j_hole_act,k_hole_act)
|
||||
|
||||
! first hole
|
||||
ispin = hole_list_practical(1,1)
|
||||
i_hole_act = hole_list_practical(2,1)
|
||||
! second hole
|
||||
jspin = hole_list_practical(1,2)
|
||||
j_hole_act = hole_list_practical(2,2)
|
||||
! third hole
|
||||
kspin = hole_list_practical(1,3)
|
||||
k_hole_act = hole_list_practical(2,3)
|
||||
delta_e_act += three_anhil(i_hole_act,j_hole_act,k_hole_act,ispin,jspin,kspin)
|
||||
|
||||
else if (n_holes_act == 0 .and. n_particles_act == 3) then
|
||||
! i_particle_act = particles_active_list_spin_traced(1)
|
||||
! j_particle_act = particles_active_list_spin_traced(2)
|
||||
! k_particle_act = particles_active_list_spin_traced(3)
|
||||
! delta_e_act += three_creat_spin_trace(i_particle_act,j_particle_act,k_particle_act)
|
||||
! first particle
|
||||
ispin = particle_list_practical(1,1)
|
||||
i_particle_act = particle_list_practical(2,1)
|
||||
! second particle
|
||||
jspin = particle_list_practical(1,2)
|
||||
j_particle_act = particle_list_practical(2,2)
|
||||
! second particle
|
||||
kspin = particle_list_practical(1,3)
|
||||
k_particle_act = particle_list_practical(2,3)
|
||||
|
||||
delta_e_act += three_creat(i_particle_act,j_particle_act,k_particle_act,ispin,jspin,kspin)
|
||||
|
||||
endif
|
||||
|
||||
!print*, 'one_anhil_spin_trace'
|
||||
!print*, one_anhil_spin_trace(1), one_anhil_spin_trace(2)
|
||||
|
||||
|
||||
delta_e_final = delta_e_act + delta_e_inactive - delta_e_virt
|
||||
!if(delta_e_final .le. -100d0.or.delta_e_final > 0.d0 .or. delta_e_final == 0.d0)then
|
||||
!if(delta_e_final == 0.d0)then
|
||||
call debug_det(det_1,N_int)
|
||||
call debug_det(det_2,N_int)
|
||||
print*, 'n_holes_act,n_particles_act'
|
||||
print*, n_holes_act,n_particles_act
|
||||
print*, 'delta_e_act,delta_e_inactive,delta_e_vir'
|
||||
print*, delta_e_act,delta_e_inactive,delta_e_virt
|
||||
delta_e_final = -1000.d0
|
||||
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
|
||||
|
||||
@ -7,12 +7,6 @@ use bitmasks
|
||||
integer :: i
|
||||
number_of_holes = 0
|
||||
|
||||
do i = 1, N_int
|
||||
number_of_holes = number_of_holes &
|
||||
+ popcnt( xor( iand(reunion_of_core_inact_bitmask(i,1), xor(key_in(i,1),iand(key_in(i,1),cas_bitmask(i,1,1)))), reunion_of_core_inact_bitmask(i,1)) )&
|
||||
+ popcnt( xor( iand(reunion_of_core_inact_bitmask(i,1), xor(key_in(i,2),iand(key_in(i,2),cas_bitmask(i,1,1)))), reunion_of_core_inact_bitmask(i,1)) )
|
||||
enddo
|
||||
return
|
||||
if(N_int == 1)then
|
||||
number_of_holes = number_of_holes &
|
||||
+ popcnt( xor( iand(reunion_of_core_inact_bitmask(1,1), xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1)))), reunion_of_core_inact_bitmask(1,1)) )&
|
||||
|
||||
@ -1,5 +1,22 @@
|
||||
BEGIN_PROVIDER [ double precision, one_body_dm_mo_alpha, (mo_tot_num_align,mo_tot_num) ]
|
||||
&BEGIN_PROVIDER [ double precision, one_body_dm_mo_beta, (mo_tot_num_align,mo_tot_num) ]
|
||||
BEGIN_PROVIDER [ double precision, one_body_dm_mo_alpha_average, (mo_tot_num_align,mo_tot_num) ]
|
||||
&BEGIN_PROVIDER [ double precision, one_body_dm_mo_beta_average, (mo_tot_num_align,mo_tot_num) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Alpha and beta one-body density matrix for each state
|
||||
END_DOC
|
||||
|
||||
integer :: i
|
||||
|
||||
one_body_dm_mo_alpha_average = 0.d0
|
||||
one_body_dm_mo_beta_average = 0.d0
|
||||
do i = 1,N_states
|
||||
one_body_dm_mo_alpha_average(:,:) += one_body_dm_mo_alpha(:,:,i) * state_average_weight(i)
|
||||
one_body_dm_mo_beta_average(:,:) += one_body_dm_mo_beta(:,:,i) * state_average_weight(i)
|
||||
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
|
||||
BEGIN_DOC
|
||||
! Alpha and beta one-body density matrix for each state
|
||||
@ -11,36 +28,31 @@
|
||||
double precision :: phase
|
||||
integer :: h1,h2,p1,p2,s1,s2, degree
|
||||
integer :: exc(0:2,2,2),n_occ(2)
|
||||
double precision, allocatable :: tmp_a(:,:), tmp_b(:,:)
|
||||
double precision, allocatable :: tmp_a(:,:,:), tmp_b(:,:,:)
|
||||
|
||||
if(only_single_double_dm)then
|
||||
print*,'ONLY DOUBLE DM'
|
||||
one_body_dm_mo_alpha = one_body_single_double_dm_mo_alpha
|
||||
one_body_dm_mo_beta = one_body_single_double_dm_mo_beta
|
||||
else
|
||||
one_body_dm_mo_alpha = 0.d0
|
||||
one_body_dm_mo_beta = 0.d0
|
||||
!$OMP PARALLEL DEFAULT(NONE) &
|
||||
!$OMP PRIVATE(j,k,l,m,occ,ck, cl, ckl,phase,h1,h2,p1,p2,s1,s2, degree,exc, &
|
||||
!$OMP tmp_a, tmp_b, n_occ)&
|
||||
!$OMP SHARED(psi_det,psi_coef,N_int,N_states,state_average_weight,elec_alpha_num,&
|
||||
!$OMP SHARED(psi_det,psi_coef,N_int,N_states,elec_alpha_num,&
|
||||
!$OMP elec_beta_num,one_body_dm_mo_alpha,one_body_dm_mo_beta,N_det,mo_tot_num_align,&
|
||||
!$OMP mo_tot_num)
|
||||
allocate(tmp_a(mo_tot_num_align,mo_tot_num), tmp_b(mo_tot_num_align,mo_tot_num) )
|
||||
allocate(tmp_a(mo_tot_num_align,mo_tot_num,N_states), tmp_b(mo_tot_num_align,mo_tot_num,N_states) )
|
||||
tmp_a = 0.d0
|
||||
tmp_b = 0.d0
|
||||
!$OMP DO SCHEDULE(dynamic)
|
||||
do k=1,N_det
|
||||
call bitstring_to_list_ab(psi_det(1,1,k), occ, n_occ, N_int)
|
||||
do m=1,N_states
|
||||
ck = psi_coef(k,m)*psi_coef(k,m) * state_average_weight(m)
|
||||
ck = psi_coef(k,m)*psi_coef(k,m)
|
||||
do l=1,elec_alpha_num
|
||||
j = occ(l,1)
|
||||
tmp_a(j,j) += ck
|
||||
tmp_a(j,j,m) += ck
|
||||
enddo
|
||||
do l=1,elec_beta_num
|
||||
j = occ(l,2)
|
||||
tmp_b(j,j) += ck
|
||||
tmp_b(j,j,m) += ck
|
||||
enddo
|
||||
enddo
|
||||
do l=1,k-1
|
||||
@ -51,28 +63,27 @@
|
||||
call get_mono_excitation(psi_det(1,1,k),psi_det(1,1,l),exc,phase,N_int)
|
||||
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
||||
do m=1,N_states
|
||||
ckl = psi_coef(k,m) * psi_coef(l,m) * phase * state_average_weight(m)
|
||||
ckl = psi_coef(k,m) * psi_coef(l,m) * phase
|
||||
if (s1==1) then
|
||||
tmp_a(h1,p1) += ckl
|
||||
tmp_a(p1,h1) += ckl
|
||||
tmp_a(h1,p1,m) += ckl
|
||||
tmp_a(p1,h1,m) += ckl
|
||||
else
|
||||
tmp_b(h1,p1) += ckl
|
||||
tmp_b(p1,h1) += ckl
|
||||
tmp_b(h1,p1,m) += ckl
|
||||
tmp_b(p1,h1,m) += ckl
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO NOWAIT
|
||||
!$OMP CRITICAL
|
||||
one_body_dm_mo_alpha = one_body_dm_mo_alpha + tmp_a
|
||||
one_body_dm_mo_alpha(:,:,:) = one_body_dm_mo_alpha(:,:,:) + tmp_a(:,:,:)
|
||||
!$OMP END CRITICAL
|
||||
!$OMP CRITICAL
|
||||
one_body_dm_mo_beta = one_body_dm_mo_beta + tmp_b
|
||||
one_body_dm_mo_beta(:,:,:) = one_body_dm_mo_beta(:,:,:) + tmp_b(:,:,:)
|
||||
!$OMP END CRITICAL
|
||||
deallocate(tmp_a,tmp_b)
|
||||
!$OMP END PARALLEL
|
||||
|
||||
endif
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, one_body_single_double_dm_mo_alpha, (mo_tot_num_align,mo_tot_num) ]
|
||||
@ -163,7 +174,7 @@ BEGIN_PROVIDER [ double precision, one_body_dm_mo, (mo_tot_num_align,mo_tot_num)
|
||||
BEGIN_DOC
|
||||
! One-body density matrix
|
||||
END_DOC
|
||||
one_body_dm_mo = one_body_dm_mo_alpha + one_body_dm_mo_beta
|
||||
one_body_dm_mo = one_body_dm_mo_alpha_average + one_body_dm_mo_beta_average
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, one_body_spin_density_mo, (mo_tot_num_align,mo_tot_num) ]
|
||||
@ -171,7 +182,7 @@ BEGIN_PROVIDER [ double precision, one_body_spin_density_mo, (mo_tot_num_align,m
|
||||
BEGIN_DOC
|
||||
! rho(alpha) - rho(beta)
|
||||
END_DOC
|
||||
one_body_spin_density_mo = one_body_dm_mo_alpha - one_body_dm_mo_beta
|
||||
one_body_spin_density_mo = one_body_dm_mo_alpha_average - one_body_dm_mo_beta_average
|
||||
END_PROVIDER
|
||||
|
||||
subroutine set_natural_mos
|
||||
@ -246,8 +257,8 @@ END_PROVIDER
|
||||
do l = 1, ao_num
|
||||
do i = 1, mo_tot_num
|
||||
do j = 1, mo_tot_num
|
||||
mo_alpha = one_body_dm_mo_alpha(j,i)
|
||||
mo_beta = one_body_dm_mo_beta(j,i)
|
||||
mo_alpha = one_body_dm_mo_alpha_average(j,i)
|
||||
mo_beta = one_body_dm_mo_beta_average(j,i)
|
||||
! 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
|
||||
|
||||
Loading…
Reference in New Issue
Block a user