LCPQ
/
quantum_package
mirror of https://github.com/LCPQ/quantum_package
10
0
Fork 0
Code Issues Releases Wiki Activity

debugging FOBOCI

This commit is contained in:
Emmanuel Giner 2017-03-16 21:21:27 +01:00
parent 4e0c71df10
commit a72b890b92
22 changed files with 514 additions and 128 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
config/gfortran.cfg
View File

@ -35,14 +35,14 @@ OPENMP : 1 ; Append OpenMP flags
# -ffast-math and the Fortran-specific
# -fno-protect-parens and -fstack-arrays.
[OPT]
FCFLAGS : -Ofast
FCFLAGS :
# Profiling flags
#################
#
[PROFILE]
FC : -p -g
FCFLAGS : -Ofast
FCFLAGS :
# Debugging flags
#################

2
plugins/FCIdump/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Determinants Davidson
Determinants Davidson core_integrals

48
plugins/FCIdump/fcidump.irp.f
View File

@ -1,21 +1,25 @@
program fcidump
implicit none
character*(128) :: output
integer :: i_unit_output,getUnitAndOpen
output=trim(ezfio_filename)//'.FCIDUMP'
i_unit_output = getUnitAndOpen(output,'w')
integer :: i,j,k,l
integer :: ii(8), jj(8), kk(8),ll(8)
integer :: i1,j1,k1,l1
integer :: i2,j2,k2,l2
integer*8 :: m
character*(2), allocatable :: A(:)
print *, '&FCI NORB=', mo_tot_num, ', NELEC=', elec_num, &
write(i_unit_output,*) '&FCI NORB=', n_act_orb, ', NELEC=', elec_num-n_core_orb*2, &
', MS2=', (elec_alpha_num-elec_beta_num), ','
allocate (A(mo_tot_num))
allocate (A(n_act_orb))
A = '1,'
print *, 'ORBSYM=', (A(i), i=1,mo_tot_num)
print *,'ISYM=0,'
print *,'/'
write(i_unit_output,*) 'ORBSYM=', (A(i), i=1,n_act_orb)
write(i_unit_output,*) 'ISYM=0,'
write(i_unit_output,*) '/'
deallocate(A)
integer*8 :: i8, k1
integer(key_kind), allocatable :: keys(:)
double precision, allocatable :: values(:)
integer(cache_map_size_kind) :: n_elements, n_elements_max
@ -23,14 +27,18 @@ program fcidump
double precision :: get_mo_bielec_integral, integral
do l=1,mo_tot_num
do k=1,mo_tot_num
do j=l,mo_tot_num
do i=k,mo_tot_num
if (i>=j) then
integral = get_mo_bielec_integral(i,j,k,l,mo_integrals_map)
do l=1,n_act_orb
l1 = list_act(l)
do k=1,n_act_orb
k1 = list_act(k)
do j=l,n_act_orb
j1 = list_act(j)
do i=k,n_act_orb
i1 = list_act(i)
if (i1>=j1) then
integral = get_mo_bielec_integral(i1,j1,k1,l1,mo_integrals_map)
if (dabs(integral) > mo_integrals_threshold) then
print *, integral, i,k,j,l
write(i_unit_output,*) integral, i,k,j,l
endif
end if
enddo
@ -38,13 +46,15 @@ program fcidump
enddo
enddo
do j=1,mo_tot_num
do i=j,mo_tot_num
integral = mo_mono_elec_integral(i,j)
do j=1,n_act_orb
j1 = list_act(j)
do i=j,n_act_orb
i1 = list_act(i)
integral = mo_mono_elec_integral(i1,j1) + core_fock_operator(i1,j1)
if (dabs(integral) > mo_integrals_threshold) then
print *, integral, i,j,0,0
write(i_unit_output,*) integral, i,j,0,0
endif
enddo
enddo
print *, 0.d0, 0, 0, 0, 0
write(i_unit_output,*) core_energy, 0, 0, 0, 0
end

2
plugins/FOBOCI/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Perturbation Selectors_no_sorted Hartree_Fock Davidson CISD
Perturbation Selectors_no_sorted SCF_density Davidson CISD

1
plugins/FOBOCI/all_singles.irp.f
View File

@ -48,6 +48,7 @@ subroutine all_single(e_pt2)
print*,'-----------------------'
print*,'i = ',i
call H_apply_just_mono(pt2, norm_pert, H_pert_diag, N_st)
call make_s2_eigenfunction_first_order
call diagonalize_CI
print*,'N_det = ',N_det
print*,'E = ',CI_energy(1)

133
plugins/FOBOCI/create_1h_or_1p.irp.f
View File

@ -29,21 +29,13 @@ subroutine create_restart_and_1h(i_hole)
enddo
enddo
enddo
integer :: N_det_old
N_det_old = N_det
N_det += n_new_det
allocate (new_det(N_int,2,n_new_det))
if (psi_det_size < N_det) then
psi_det_size = N_det
TOUCH psi_det_size
endif
do i = 1, N_det_old
do k = 1, N_int
psi_det(k,1,i) = old_psi_det(k,1,i)
psi_det(k,2,i) = old_psi_det(k,2,i)
enddo
enddo
logical, allocatable :: duplicate(:)
allocate (new_det(N_int,2,n_new_det),duplicate(n_new_det))
n_new_det = 0
do j = 1, n_act_orb
@ -58,19 +50,56 @@ subroutine create_restart_and_1h(i_hole)
if(i_ok .ne. 1)cycle
n_new_det +=1
do k = 1, N_int
psi_det(k,1,n_det_old+n_new_det) = key_tmp(k,1)
psi_det(k,2,n_det_old+n_new_det) = key_tmp(k,2)
new_det(k,1,n_new_det) = key_tmp(k,1)
new_det(k,2,n_new_det) = key_tmp(k,2)
enddo
psi_coef(n_det_old+n_new_det,:) = 0.d0
enddo
enddo
enddo
SOFT_TOUCH N_det psi_det psi_coef
logical :: found_duplicates
if(n_act_orb.gt.1)then
call remove_duplicates_in_psi_det(found_duplicates)
integer :: i_test
duplicate = .False.
do i = 1, n_new_det
if(duplicate(i))cycle
do j = i+1, n_new_det
i_test = 0
do ispin =1 ,2
do k = 1, N_int
i_test += popcnt(xor(new_det(k,ispin,i),new_det(k,ispin,j)))
enddo
enddo
if(i_test.eq.0)then
duplicate(j) = .True.
endif
enddo
enddo
integer :: n_new_det_unique
n_new_det_unique = 0
print*, 'uniq det'
do i = 1, n_new_det
if(.not.duplicate(i))then
n_new_det_unique += 1
endif
enddo
print*, n_new_det_unique
N_det += n_new_det_unique
if (psi_det_size < N_det) then
psi_det_size = N_det
TOUCH psi_det_size
endif
do i = 1, n_new_det_unique
do ispin = 1, 2
do k = 1, N_int
psi_det(k,ispin,N_det_old+i) = new_det(k,ispin,i)
enddo
enddo
psi_coef(N_det_old+i,:) = 0.d0
enddo
SOFT_TOUCH N_det psi_det psi_coef
deallocate (new_det,duplicate)
end
subroutine create_restart_and_1p(i_particle)
@ -107,18 +136,8 @@ subroutine create_restart_and_1p(i_particle)
integer :: N_det_old
N_det_old = N_det
N_det += n_new_det
allocate (new_det(N_int,2,n_new_det))
if (psi_det_size < N_det) then
psi_det_size = N_det
TOUCH psi_det_size
endif
do i = 1, N_det_old
do k = 1, N_int
psi_det(k,1,i) = old_psi_det(k,1,i)
psi_det(k,2,i) = old_psi_det(k,2,i)
enddo
enddo
logical, allocatable :: duplicate(:)
allocate (new_det(N_int,2,n_new_det),duplicate(n_new_det))
n_new_det = 0
do j = 1, n_act_orb
@ -133,17 +152,59 @@ subroutine create_restart_and_1p(i_particle)
if(i_ok .ne. 1)cycle
n_new_det +=1
do k = 1, N_int
psi_det(k,1,n_det_old+n_new_det) = key_tmp(k,1)
psi_det(k,2,n_det_old+n_new_det) = key_tmp(k,2)
new_det(k,1,n_new_det) = key_tmp(k,1)
new_Det(k,2,n_new_det) = key_tmp(k,2)
enddo
psi_coef(n_det_old+n_new_det,:) = 0.d0
enddo
enddo
enddo
integer :: i_test
duplicate = .False.
do i = 1, n_new_det
if(duplicate(i))cycle
call debug_det(new_det(1,1,i),N_int)
do j = i+1, n_new_det
i_test = 0
call debug_det(new_det(1,1,j),N_int)
do ispin =1 ,2
do k = 1, N_int
i_test += popcnt(xor(new_det(k,ispin,i),new_det(k,ispin,j)))
enddo
enddo
if(i_test.eq.0)then
duplicate(j) = .True.
endif
enddo
enddo
integer :: n_new_det_unique
n_new_det_unique = 0
print*, 'uniq det'
do i = 1, n_new_det
if(.not.duplicate(i))then
n_new_det_unique += 1
endif
enddo
print*, n_new_det_unique
N_det += n_new_det_unique
if (psi_det_size < N_det) then
psi_det_size = N_det
TOUCH psi_det_size
endif
do i = 1, n_new_det_unique
do ispin = 1, 2
do k = 1, N_int
psi_det(k,ispin,N_det_old+i) = new_det(k,ispin,i)
enddo
enddo
psi_coef(N_det_old+i,:) = 0.d0
enddo
SOFT_TOUCH N_det psi_det psi_coef
logical :: found_duplicates
call remove_duplicates_in_psi_det(found_duplicates)
deallocate (new_det,duplicate)
end
subroutine create_restart_1h_1p(i_hole,i_part)

5
plugins/FOBOCI/density_matrix.irp.f
View File

@ -32,6 +32,11 @@
psi_coef_generators_restart(i,1) = psi_coef_generators_restart(i,1) * inv_coef_ref_generators_restart
norm_generators_restart += psi_coef_generators_restart(i,1)**2
enddo
double precision :: inv_norm
inv_norm = 1.d0/dsqrt(norm_generators_restart)
do i = 1, N_det_generators_restart
psi_coef_generators_restart(i,1) = psi_coef_generators_restart(i,1) * inv_norm
enddo
one_body_dm_mo_alpha_generators_restart = 0.d0

7
plugins/FOBOCI/dress_simple.irp.f
View File

@ -175,6 +175,10 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener
endif
do j = 1, Ndet_generators
call i_h_j(psi_det_generators_input(1,1,j),psi_det_generators_input(1,1,i),N_int,hij) ! Fill the zeroth order H matrix
if(i==j)then
call debug_det(psi_det_generators_input(1,1,i),N_int)
print*, hij
endif
dressed_H_matrix(i,j) = hij
enddo
enddo
@ -234,6 +238,7 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener
do i = 1, N_states
i_state(i) = i
E_ref(i) = eigvalues(i)
print*, 'E_ref(i)',E_ref(i)
enddo
endif
do i = 1,N_states
@ -287,7 +292,7 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener
accu += eigvectors(j,i) * psi_coef_ref(j,k)
enddo
print*,'accu = ',accu
if(dabs(accu).ge.0.72d0)then
if(dabs(accu).ge.0.60d0)then
i_good_state(0) +=1
i_good_state(i_good_state(0)) = i
endif

6
plugins/FOBOCI/fobo_scf.irp.f
View File

@ -15,8 +15,6 @@ end
subroutine run_prepare
implicit none
! no_oa_or_av_opt = .False.
! touch no_oa_or_av_opt
call damping_SCF
call diag_inactive_virt_and_update_mos
end
@ -28,7 +26,7 @@ subroutine routine_fobo_scf
print*,''
character*(64) :: label
label = "Natural"
do i = 1, 5
do i = 1, 1
print*,'*******************************************************************************'
print*,'*******************************************************************************'
print*,'FOBO-SCF Iteration ',i
@ -54,7 +52,7 @@ subroutine routine_fobo_scf
endif
call FOBOCI_lmct_mlct_old_thr(i)
call save_osoci_natural_mos
call damping_SCF
! call damping_SCF
call diag_inactive_virt_and_update_mos
call clear_mo_map
call provide_properties

7
plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f
View File

@ -55,6 +55,10 @@ subroutine FOBOCI_lmct_mlct_old_thr(iter)
call create_restart_and_1h(i_hole_osoci)
call set_generators_to_psi_det
print*,'Passed set generators'
integer :: m
do m = 1, N_det_generators
call debug_det(psi_det_generators(1,1,m),N_int)
enddo
call set_bitmask_particl_as_input(reunion_of_bitmask)
call set_bitmask_hole_as_input(reunion_of_bitmask)
double precision :: e_pt2
@ -82,7 +86,7 @@ subroutine FOBOCI_lmct_mlct_old_thr(iter)
call set_bitmask_particl_as_input(reunion_of_bitmask)
call set_bitmask_hole_as_input(reunion_of_bitmask)
call all_single(e_pt2)
call make_s2_eigenfunction_first_order
! call make_s2_eigenfunction_first_order
threshold_davidson = 1.d-6
soft_touch threshold_davidson davidson_criterion
call diagonalize_ci
@ -541,7 +545,6 @@ subroutine FOBOCI_lmct_mlct_old_thr_restart(iter)
call print_generators_bitmasks_holes
! Impose that only the active part can be reached
call set_bitmask_hole_as_input(unpaired_bitmask)
!!! call all_single_h_core
call create_restart_and_1p(i_particl_osoci)
!!! ! Update the generators
call set_generators_to_psi_det

110
plugins/FOBOCI/routines_foboci.irp.f
View File

@ -13,6 +13,8 @@ subroutine set_intermediate_normalization_lmct_old(norm,i_hole)
integer :: n_good_hole
logical,allocatable :: is_a_ref_det(:)
allocate(index_one_hole(n_det),index_one_hole_one_p(n_det),index_two_hole_one_p(N_det),index_two_hole(N_det),index_one_p(N_det),is_a_ref_det(N_det))
double precision, allocatable :: local_norm(:)
allocate(local_norm(N_states))
n_one_hole = 0
n_one_hole_one_p = 0
@ -30,7 +32,6 @@ subroutine set_intermediate_normalization_lmct_old(norm,i_hole)
do k = 1, N_states
inv_coef_ref_generators_restart(k) = 1.d0/psi_coef(i,k)
enddo
! cycle
endif
! Find all the determinants present in the reference wave function
@ -59,10 +60,8 @@ subroutine set_intermediate_normalization_lmct_old(norm,i_hole)
enddo
endif
enddo
!do k = 1, N_det
! call debug_det(psi_det(1,1,k),N_int)
! print*,'k,coef = ',k,psi_coef(k,1)/psi_coef(index_ref_generators_restart,1)
!enddo
print*,''
print*,'n_good_hole = ',n_good_hole
do k = 1,N_states
@ -72,27 +71,37 @@ subroutine set_intermediate_normalization_lmct_old(norm,i_hole)
enddo
print*,''
enddo
norm = 0.d0
! Set the wave function to the intermediate normalization
! Set the wave function to the intermediate normalization
do k = 1, N_states
do i = 1, N_det
psi_coef(i,k) = psi_coef(i,k) * inv_coef_ref_generators_restart(k)
enddo
enddo
norm = 0.d0
do k = 1,N_states
print*,'state ',k
do i = 1, N_det
!! print*,'psi_coef(i_ref) = ',psi_coef(i,1)
if (is_a_ref_det(i))then
print*,'i,psi_coef_ref = ',psi_coef(i,k)
cycle
endif
norm(k) += psi_coef(i,k) * psi_coef(i,k)
enddo
print*,'norm = ',norm(k)
enddo
do k =1, N_states
local_norm(k) = 1.d0 / dsqrt(norm(k))
enddo
do k = 1,N_states
do i = 1, N_det
psi_coef(i,k) = psi_coef(i,k) * local_norm(k)
enddo
enddo
deallocate(index_one_hole,index_one_hole_one_p,index_two_hole_one_p,index_two_hole,index_one_p,is_a_ref_det)
deallocate(local_norm)
soft_touch psi_coef
end
@ -117,6 +126,8 @@ subroutine set_intermediate_normalization_mlct_old(norm,i_particl)
integer :: i_count
allocate(index_one_hole(n_det),index_one_hole_one_p(n_det),index_two_hole_one_p(N_det),index_two_hole(N_det),index_one_p(N_det),is_a_ref_det(N_det))
allocate(index_one_hole_two_p(n_det))
double precision, allocatable :: local_norm(:)
allocate(local_norm(N_states))
n_one_hole = 0
n_one_hole_one_p = 0
@ -185,20 +196,29 @@ subroutine set_intermediate_normalization_mlct_old(norm,i_particl)
psi_coef(i,k) = psi_coef(i,k) * inv_coef_ref_generators_restart(k)
enddo
enddo
do k = 1, N_states
norm = 0.d0
do k = 1,N_states
print*,'state ',k
do i = 1, N_det
!! print*,'i = ',i, psi_coef(i,1)
if (is_a_ref_det(i))then
print*,'i,psi_coef_ref = ',psi_coef(i,k)
cycle
endif
norm(k) += psi_coef(i,k) * psi_coef(i,k)
enddo
print*,'norm = ',norm
print*,'norm = ',norm(k)
enddo
do k =1, N_states
local_norm(k) = 1.d0 / dsqrt(norm(k))
enddo
do k = 1,N_states
do i = 1, N_det
psi_coef(i,k) = psi_coef(i,k) * local_norm(k)
enddo
enddo
soft_touch psi_coef
deallocate(index_one_hole,index_one_hole_one_p,index_two_hole_one_p,index_two_hole,index_one_p,is_a_ref_det)
deallocate(local_norm)
end
@ -210,12 +230,60 @@ subroutine update_density_matrix_osoci
END_DOC
integer :: i,j
integer :: iorb,jorb
! active <--> inactive block
do i = 1, mo_tot_num
do j = 1, mo_tot_num
one_body_dm_mo_alpha_osoci(i,j) = one_body_dm_mo_alpha_osoci(i,j) + (one_body_dm_mo_alpha_average(i,j) - one_body_dm_mo_alpha_generators_restart(i,j))
one_body_dm_mo_beta_osoci(i,j) = one_body_dm_mo_beta_osoci(i,j) + (one_body_dm_mo_beta_average(i,j) - one_body_dm_mo_beta_generators_restart(i,j))
one_body_dm_mo_alpha_osoci(i,j) += one_body_dm_mo_alpha_average(i,j) - one_body_dm_mo_alpha_generators_restart(i,j)
one_body_dm_mo_beta_osoci(i,j) += one_body_dm_mo_beta_average(i,j) - one_body_dm_mo_beta_generators_restart(i,j)
enddo
enddo
!do i = 1, n_act_orb
! iorb = list_act(i)
! do j = 1, n_inact_orb
! jorb = list_inact(j)
! one_body_dm_mo_alpha_osoci(iorb,jorb)+= one_body_dm_mo_alpha_average(iorb,jorb)
! one_body_dm_mo_alpha_osoci(jorb,iorb)+= one_body_dm_mo_alpha_average(jorb,iorb)
! one_body_dm_mo_beta_osoci(iorb,jorb) += one_body_dm_mo_beta_average(iorb,jorb)
! one_body_dm_mo_beta_osoci(jorb,iorb) += one_body_dm_mo_beta_average(jorb,iorb)
! enddo
!enddo
!! active <--> virt block
!do i = 1, n_act_orb
! iorb = list_act(i)
! do j = 1, n_virt_orb
! jorb = list_virt(j)
! one_body_dm_mo_alpha_osoci(iorb,jorb)+= one_body_dm_mo_alpha_average(iorb,jorb)
! one_body_dm_mo_alpha_osoci(jorb,iorb)+= one_body_dm_mo_alpha_average(jorb,iorb)
! one_body_dm_mo_beta_osoci(iorb,jorb) += one_body_dm_mo_beta_average(iorb,jorb)
! one_body_dm_mo_beta_osoci(jorb,iorb) += one_body_dm_mo_beta_average(jorb,iorb)
! enddo
!enddo
!! virt <--> virt block
!do j = 1, n_virt_orb
! jorb = list_virt(j)
! one_body_dm_mo_alpha_osoci(jorb,jorb)+= one_body_dm_mo_alpha_average(jorb,jorb)
! one_body_dm_mo_beta_osoci(jorb,jorb) += one_body_dm_mo_beta_average(jorb,jorb)
!enddo
!! inact <--> inact block
!do j = 1, n_inact_orb
! jorb = list_inact(j)
! one_body_dm_mo_alpha_osoci(jorb,jorb) -= one_body_dm_mo_alpha_average(jorb,jorb)
! one_body_dm_mo_beta_osoci(jorb,jorb) -= one_body_dm_mo_beta_average(jorb,jorb)
!enddo
double precision :: accu_alpha, accu_beta
accu_alpha = 0.d0
accu_beta = 0.d0
do i = 1, mo_tot_num
accu_alpha += one_body_dm_mo_alpha_osoci(i,i)
accu_beta += one_body_dm_mo_beta_osoci(i,i)
! write(*,'(I3,X,100(F16.10,X))') i,one_body_dm_mo_alpha_osoci(i,i),one_body_dm_mo_beta_osoci(i,i),one_body_dm_mo_alpha_osoci(i,i)+one_body_dm_mo_beta_osoci(i,i)
enddo
print*, 'accu_alpha/beta',accu_alpha,accu_beta
end
@ -263,6 +331,12 @@ subroutine initialize_density_matrix_osoci
implicit none
one_body_dm_mo_alpha_osoci = one_body_dm_mo_alpha_generators_restart
one_body_dm_mo_beta_osoci = one_body_dm_mo_beta_generators_restart
integer :: i
print*, '8*********************'
print*, 'initialize_density_matrix_osoci'
do i = 1, mo_tot_num
print*,one_body_dm_mo_alpha_osoci(i,i),one_body_dm_mo_alpha_generators_restart(i,i)
enddo
end
subroutine rescale_density_matrix_osoci(norm)
@ -438,6 +512,10 @@ subroutine save_osoci_natural_mos
endif
enddo
enddo
print*, 'test'
print*, 'test'
print*, 'test'
print*, 'test'
do i = 1, mo_tot_num
do j = i+1, mo_tot_num
if(dabs(tmp(i,j)).le.threshold_fobo_dm)then
@ -445,7 +523,9 @@ subroutine save_osoci_natural_mos
tmp(j,i) = 0.d0
endif
enddo
print*, tmp(i,i)
enddo
label = "Natural"

10
plugins/MRPT/MRPT_Utils.main.irp.f
View File

@ -18,11 +18,11 @@ subroutine routine_3
print *, 'N_states = ', N_states
do i = 1, N_States
print*,'State',i
write(*,'(A12,X,I3,A3,XX,F16.10)') ' PT2 ', i,' = ', second_order_pt_new(i)
write(*,'(A12,X,I3,A3,XX,F16.09)') ' E ', i,' = ', psi_ref_average_value(i)
write(*,'(A12,X,I3,A3,XX,F16.09)') ' E+PT2 ', i,' = ', psi_ref_average_value(i)+second_order_pt_new(i)
write(*,'(A12,X,I3,A3,XX,F16.09)') ' E dressed ', i,' = ', CI_dressed_pt2_new_energy(i)
write(*,'(A12,X,I3,A3,XX,F16.09)') ' S^2 ', i,' = ', CI_dressed_pt2_new_eigenvectors_s2(i)
write(*,'(A12,X,I3,A3,XX,F20.16)') ' PT2 ', i,' = ', second_order_pt_new(i)
write(*,'(A12,X,I3,A3,XX,F22.16)') ' E ', i,' = ', psi_ref_average_value(i)
write(*,'(A12,X,I3,A3,XX,F22.16)') ' E+PT2 ', i,' = ', psi_ref_average_value(i)+second_order_pt_new(i)
write(*,'(A12,X,I3,A3,XX,F22.16)') ' E dressed ', i,' = ', CI_dressed_pt2_new_energy(i)
write(*,'(A12,X,I3,A3,XX,F20.16)') ' S^2 ', i,' = ', CI_dressed_pt2_new_eigenvectors_s2(i)
print*,'coef before and after'
do j = 1, N_det_ref
print*,psi_ref_coef(j,i),CI_dressed_pt2_new_eigenvectors(j,i)

8
plugins/MRPT/print_1h2p.irp.f
View File

@ -5,6 +5,12 @@ program print_1h2p
call routine
end
subroutine routine
implicit none
provide one_anhil_one_creat_inact_virt
end
subroutine routine_2
implicit none
integer :: i,j,degree
@ -27,7 +33,7 @@ subroutine routine_2
end
subroutine routine
subroutine routine_3
implicit none
double precision,allocatable :: matrix_1h2p(:,:,:)
double precision :: accu(2)

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

@ -22,6 +22,40 @@ BEGIN_PROVIDER [ double precision, energy_cas_dyall_no_exchange, (N_states)]
END_PROVIDER
BEGIN_PROVIDER [ double precision, energy_cas_dyall_no_exchange_bis, (N_states)]
implicit none
integer :: i,j
double precision :: energies(N_states)
integer(bit_kind), allocatable :: psi_in_ref(:,:,:)
allocate (psi_in_ref(N_int,2,n_det_ref))
integer(bit_kind), allocatable :: psi_in_active(:,:,:)
allocate (psi_in_active(N_int,2,n_det_ref))
double precision, allocatable :: psi_ref_coef_in(:, :)
allocate(psi_ref_coef_in(N_det_ref, N_states))
do i = 1, N_det_ref
do j = 1, N_int
psi_in_ref(j,1,i) = psi_ref(j,1,i)
psi_in_ref(j,2,i) = psi_ref(j,2,i)
psi_in_active(j,1,i) = psi_active(j,1,i)
psi_in_active(j,2,i) = psi_active(j,2,i)
enddo
do j = 1, N_states
psi_ref_coef_in(i,j) = psi_ref_coef(i,j)
enddo
enddo
do i = 1, N_states
call u0_H_dyall_u0_no_exchange_bis(energies,psi_in_ref,psi_ref_coef_in,n_det_ref,n_det_ref,n_det_ref,N_states,i)
energy_cas_dyall_no_exchange_bis(i) = energies(i)
print*, 'energy_cas_dyall(i)_no_exchange_bis', energy_cas_dyall_no_exchange_bis(i)
enddo
deallocate (psi_in_ref)
deallocate (psi_in_active)
deallocate(psi_ref_coef_in)
END_PROVIDER
BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2,N_states)]
implicit none
@ -604,6 +638,8 @@ END_PROVIDER
double precision, allocatable :: psi_in_out_coef(:,:)
use bitmasks
allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states))
integer(bit_kind), allocatable :: psi_in_active(:,:,:)
allocate (psi_in_active(N_int,2,n_det_ref))
integer :: iorb,jorb,i_ok
integer :: state_target
@ -614,6 +650,9 @@ END_PROVIDER
double precision :: thresh_norm
integer :: other_spin(2)
other_spin(1) = 2
other_spin(2) = 1
thresh_norm = 1.d-20
@ -644,14 +683,14 @@ END_PROVIDER
print*, 'pb, i_ok ne 0 !!!'
endif
call i_H_j(psi_in_out(1,1,i),psi_ref(1,1,i),N_int,hij)
integer :: exc(0:2,2,2)
double precision :: phase
call get_mono_excitation(psi_in_out(1,1,i),psi_ref(1,1,i),exc,phase,N_int)
do j = 1, n_states
double precision :: coef,contrib
coef = psi_ref_coef(i,j) !* psi_ref_coef(i,j)
psi_in_out_coef(i,j) = coef * hij
if(orb_i == 5 .and. orb_v == 20)then
! if(orb_i == 2 .and. orb_v == 6 )then
psi_in_out_coef(i,j) = psi_ref_coef(i,j)* hij * phase
! if(orb_i == 5 .and. orb_v == 20)then
if(orb_i == 2 .and. orb_v == 6 )then
print*, i, ispin
print*, coef * hij,coef,hij
endif
norm(j,ispin) += psi_in_out_coef(i,j) * psi_in_out_coef(i,j)
enddo
@ -663,22 +702,31 @@ END_PROVIDER
one_anhil_one_creat_inact_virt_norm(iorb,vorb,j,ispin) = 0.d0
else
norm_no_inv(j,ispin) = norm(j,ispin)
one_anhil_one_creat_inact_virt_norm(iorb,vorb,j,ispin) = 1.d0 / norm(j,ispin)
! one_anhil_one_creat_inact_virt_norm(iorb,vorb,j,ispin) = 1.d0 / norm(j,ispin)
norm(j,ispin) = 1.d0/dsqrt(norm(j,ispin))
if(orb_i == 5 .and. orb_v == 20)then
! if(orb_i == 2 .and. orb_v == 6 )then
! if(orb_i == 5 .and. orb_v == 20)then
if(orb_i == 2 .and. orb_v == 6 )then
print*,ispin ,norm(j,ispin)
endif
endif
enddo
integer :: iorb_annil,hole_particle,spin_exc,orb
double precision :: norm_out_bis(N_states)
do i = 1, N_det_ref
do j = 1, N_states
psi_in_out_coef(i,j) = psi_in_out_coef(i,j) * norm(j,ispin)
norm_bis(j,ispin) += psi_in_out_coef(i,j) * psi_in_out_coef(i,j)
enddo
enddo
do i = 1, N_det_ref
do j = 1, N_int
! psi_in_out(j,1,i) = iand(psi_in_out(j,1,i),cas_bitmask(j,1,1))
! psi_in_out(j,2,i) = iand(psi_in_out(j,2,i),cas_bitmask(j,1,1))
psi_in_out(j,1,i) = psi_active(j,1,i)
psi_in_out(j,2,i) = psi_active(j,2,i)
! psi_in_out(j,1,i) = psi_ref(j,1,i)
! psi_in_out(j,2,i) = psi_ref(j,2,i)
enddo
enddo
do state_target = 1, N_states
@ -686,29 +734,35 @@ END_PROVIDER
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_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
energies_alpha_beta(state_target, ispin) += energies(state_target)
if(orb_i == 5 .and. orb_v == 20)then
! if(orb_i == 2 .and. orb_v == 6 )then
print*, ispin, energy_cas_dyall_no_exchange(1) , energies_alpha_beta(state_target, ispin)
print*, ispin, energy_cas_dyall_no_exchange(1) - energies_alpha_beta(state_target, ispin)
endif
endif
enddo
enddo ! ispin
do state_target = 1, N_states
if((norm_no_inv(state_target,1) + norm_no_inv(state_target,2)) .ne. 0.d0)then
one_anhil_one_creat_inact_virt(iorb,vorb,state_target) = energy_cas_dyall_no_exchange(state_target) - &
( energies_alpha_beta(state_target,1) + energies_alpha_beta(state_target,2) ) &
! 0.5d0 * (energies_alpha_beta(state_target,1) + energies_alpha_beta(state_target,2))
( energies_alpha_beta(state_target,1) + energies_alpha_beta(state_target,2) ) &
/( norm_bis(state_target,1) + norm_bis(state_target,2) )
else
one_anhil_one_creat_inact_virt(iorb,vorb,state_target) = 0.d0
endif
if(dabs(dabs(one_anhil_one_creat_inact_virt(iorb,vorb,state_target)) - 1.30584271462d0) < 1.d-11)then
! if(dabs(dabs(one_anhil_one_creat_inact_virt(iorb,vorb,state_target)) - 1.29269686324d0) < 1.d-11)then
print*, ''
print*, orb_i,orb_v
print*, energy_cas_dyall_no_exchange(1) - energies_alpha_beta(state_target,1) / norm_bis(state_target,1)
print*, energy_cas_dyall_no_exchange(1) - energies_alpha_beta(state_target,2) / norm_bis(state_target,2)
print*, energy_cas_dyall_no_exchange(1) - energies_alpha_beta(state_target,1) !/ norm_bis(state_target,1)
print*, energy_cas_dyall_no_exchange(1) - energies_alpha_beta(state_target,2) !/ norm_bis(state_target,2)
print*, fock_core_inactive_total_spin_trace(orb_i,1)
print*, fock_virt_total_spin_trace(orb_v,1)
print*, one_anhil_one_creat_inact_virt(iorb,vorb,state_target)
print*, ''
endif
if(dabs(one_anhil_one_creat_inact_virt(iorb,vorb,state_target)) .gt. 1.d-10)then
write(*,'(F11.8)'), one_anhil_one_creat_inact_virt(iorb,vorb,state_target)
! if(dabs(one_anhil_one_creat_inact_virt(iorb,vorb,state_target)) .lt. 1.d-2)then
! one_anhil_one_creat_inact_virt(iorb,vorb,state_target) = 0.d0
! print*, orb_i,orb_v
! endif
endif
enddo
enddo
@ -852,9 +906,6 @@ BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_State
norm_bis = 0.d0
do ispin = 1,2
do i = 1, n_det_ref
! if(orb_a == 6 .and. orb_v == 12)then
! print*, 'i ref = ',i
! endif
do j = 1, N_int
psi_in_out(j,1,i) = psi_ref(j,1,i)
psi_in_out(j,2,i) = psi_ref(j,2,i)
@ -866,11 +917,6 @@ BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_State
enddo
else
call i_H_j(psi_in_out(1,1,i),psi_ref(1,1,i),N_int,hij)
if(orb_a == 6 .and. orb_v == 12)then
call debug_det(psi_ref(1,1,i),N_int)
call debug_det(psi_in_out(1,1,i),N_int)
print*, hij
endif
do j = 1, n_states
double precision :: contrib
psi_in_out_coef(i,j) = psi_ref_coef(i,j) * hij
@ -907,7 +953,6 @@ 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_ref,n_det_ref,n_det_ref,N_states,state_target)
call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
energies_alpha_beta(state_target, ispin) += energies(state_target)
endif
@ -915,7 +960,6 @@ BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_State
enddo ! ispin
do state_target = 1, N_states
if((norm_no_inv(state_target,1) + norm_no_inv(state_target,2)) .ne. 0.d0)then
! one_creat_virt(aorb,vorb,state_target) = energy_cas_dyall(state_target) - &
one_creat_virt(aorb,vorb,state_target) = energy_cas_dyall_no_exchange(state_target) - &
( energies_alpha_beta(state_target,1) + energies_alpha_beta(state_target,2) ) &
/( norm_bis(state_target,1) + norm_bis(state_target,2) )

55
plugins/MRPT_Utils/excitations_cas.irp.f
View File

@ -25,6 +25,7 @@ subroutine apply_exc_to_psi(orb,hole_particle,spin_exc, &
integer :: det_tmp(N_int), det_tmp_bis(N_int)
double precision :: phase
double precision :: norm_factor
! print*, orb,hole_particle,spin_exc
elec_num_tab_local = 0
do i = 1, ndet
@ -36,6 +37,7 @@ subroutine apply_exc_to_psi(orb,hole_particle,spin_exc, &
exit
endif
enddo
! print*, elec_num_tab_local(1),elec_num_tab_local(2)
if(hole_particle == 1)then
do i = 1, ndet
call set_bit_to_integer(orb,psi_in_out(1,spin_exc,i),N_int)
@ -675,6 +677,7 @@ subroutine i_H_j_dyall_no_exchange(key_i,key_j,Nint,hij)
integer :: n_occ_ab(2)
logical :: has_mipi(Nint*bit_kind_size)
double precision :: mipi(Nint*bit_kind_size)
double precision :: diag_H_mat_elem
PROVIDE mo_bielec_integrals_in_map mo_integrals_map
ASSERT (Nint > 0)
@ -771,9 +774,11 @@ subroutine i_H_j_dyall_no_exchange(key_i,key_j,Nint,hij)
endif
hij = phase*(hij + mo_mono_elec_integral(m,p) + fock_operator_active_from_core_inact(m,p) )
! hij = phase*(hij + mo_mono_elec_integral(m,p) )
case (0)
hij = diag_H_mat_elem_no_elec_check_no_exchange(key_i,Nint)
! hij = diag_H_mat_elem(key_i,Nint)
! hij = 0.d0
end select
end
@ -799,7 +804,7 @@ double precision function diag_H_mat_elem_no_elec_check_no_exchange(det_in,Nint)
! alpha - alpha
do i = 1, elec_num_tab_local(1)
iorb = occ(i,1)
diag_H_mat_elem_no_elec_check_no_exchange += mo_mono_elec_integral(iorb,iorb) !+ fock_operator_active_from_core_inact(iorb,iorb)
diag_H_mat_elem_no_elec_check_no_exchange += mo_mono_elec_integral(iorb,iorb) !+ fock_operator_active_from_core_inact(iorb,iorb)
do j = i+1, elec_num_tab_local(1)
jorb = occ(j,1)
diag_H_mat_elem_no_elec_check_no_exchange += mo_bielec_integral_jj(jorb,iorb)
@ -809,7 +814,7 @@ double precision function diag_H_mat_elem_no_elec_check_no_exchange(det_in,Nint)
! beta - beta
do i = 1, elec_num_tab_local(2)
iorb = occ(i,2)
diag_H_mat_elem_no_elec_check_no_exchange += mo_mono_elec_integral(iorb,iorb) !+ fock_operator_active_from_core_inact(iorb,iorb)
diag_H_mat_elem_no_elec_check_no_exchange += mo_mono_elec_integral(iorb,iorb) !+ fock_operator_active_from_core_inact(iorb,iorb)
do j = i+1, elec_num_tab_local(2)
jorb = occ(j,2)
diag_H_mat_elem_no_elec_check_no_exchange += mo_bielec_integral_jj(jorb,iorb)
@ -835,7 +840,7 @@ double precision function diag_H_mat_elem_no_elec_check_no_exchange(det_in,Nint)
do j = 1, n_core_inact_orb
jorb = list_core_inact(j)
diag_H_mat_elem_no_elec_check_no_exchange += 2.d0 * mo_bielec_integral_jj(jorb,iorb)
core_act_exchange(1) += - mo_bielec_integral_jj_exchange(jorb,iorb)
! core_act_exchange(1) += - mo_bielec_integral_jj_exchange(jorb,iorb)
! diag_H_mat_elem_no_elec_check_no_exchange += core_act_exchange(1)
enddo
enddo
@ -846,7 +851,7 @@ double precision function diag_H_mat_elem_no_elec_check_no_exchange(det_in,Nint)
do j = 1, n_core_inact_orb
jorb = list_core_inact(j)
diag_H_mat_elem_no_elec_check_no_exchange += 2.d0 * mo_bielec_integral_jj(jorb,iorb)
core_act_exchange(2) += - mo_bielec_integral_jj_exchange(jorb,iorb)
! core_act_exchange(2) += - mo_bielec_integral_jj_exchange(jorb,iorb)
! diag_H_mat_elem_no_elec_check_no_exchange += core_act_exchange(2)
enddo
enddo
@ -884,3 +889,45 @@ subroutine u0_H_dyall_u0_no_exchange(energies,psi_in,psi_in_coef,ndet,dim_psi_in
energies(state_target) = accu
deallocate(psi_coef_tmp)
end
!subroutine u0_H_dyall_u0_no_exchange_bis(energies,psi_in,psi_in_active,psi_in_coef,ndet,dim_psi_in,dim_psi_coef,N_states_in,state_target)
subroutine u0_H_dyall_u0_no_exchange_bis(energies,psi_in,psi_in_coef,ndet,dim_psi_in,dim_psi_coef,N_states_in,state_target)
use bitmasks
implicit none
integer, intent(in) :: N_states_in,ndet,dim_psi_in,dim_psi_coef,state_target
!integer(bit_kind), intent(in) :: psi_in(N_int,2,dim_psi_in),psi_in_active(N_int,2,dim_psi_in)
integer(bit_kind), intent(in) :: psi_in(N_int,2,dim_psi_in)
double precision, intent(in) :: psi_in_coef(dim_psi_coef,N_states_in)
double precision, intent(out) :: energies(N_states_in)
integer :: i,j
double precision :: hij,accu
energies = 0.d0
accu = 0.d0
double precision, allocatable :: psi_coef_tmp(:)
allocate(psi_coef_tmp(ndet))
do i = 1, ndet
psi_coef_tmp(i) = psi_in_coef(i,state_target)
enddo
double precision :: hij_bis,diag_H_mat_elem
do i = 1, ndet
if(psi_coef_tmp(i)==0.d0)cycle
do j = i+1, ndet
if(psi_coef_tmp(j)==0.d0)cycle
! call i_H_j_dyall_no_exchange(psi_in_active(1,1,i),psi_in_active(1,1,j),N_int,hij)
call i_H_j(psi_in(1,1,i),psi_in(1,1,j),N_int,hij)
accu += 2.d0 * psi_coef_tmp(i) * psi_coef_tmp(j) * hij
enddo
enddo
do i = 1, N_det
if(psi_coef_tmp(i)==0.d0)cycle
accu += psi_coef_tmp(i) * psi_coef_tmp(i) * diag_H_mat_elem(psi_in(1,1,i),N_int)
enddo
energies(state_target) = accu
deallocate(psi_coef_tmp)
end

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

@ -197,7 +197,7 @@
k_inact_core_orb = list_core_inact(k)
coulomb = get_mo_bielec_integral(k_inact_core_orb,iorb,k_inact_core_orb,jorb,mo_integrals_map)
exchange = get_mo_bielec_integral(k_inact_core_orb,jorb,iorb,k_inact_core_orb,mo_integrals_map)
accu += 2.d0 * coulomb - exchange
accu += 2.d0 * coulomb - exchange
enddo
fock_operator_active_from_core_inact(iorb,jorb) = accu
enddo

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

@ -122,7 +122,7 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip
enddo
else
call get_delta_e_dyall(psi_ref(1,1,index_i),tq(1,1,i_alpha),coef_array,hialpha,delta_e)
! !!!!!!!!!!!!! SHIFTED BK
!!!!!!!!!!!!! SHIFTED BK
! double precision :: hjj
! call i_h_j(tq(1,1,i_alpha),tq(1,1,i_alpha),Nint,hjj)
! delta_e(1) = CI_electronic_energy(1) - hjj
@ -141,7 +141,11 @@ 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)
!!!!!!!!!!!!!!!!!! WARNING TEST
if(index_j .ne. index_i)cycle
!!!!!!!!!!!!!!!!!! WARNING TEST
! if(index_j .ne. index_i)cycle
!!!!!!!!!!!!!!!!!! WARNING TEST
!!!!!!!!!!!!!!!!!! WARNING TEST
!!!!!!!!!!!!!!!!!! WARNING TEST
do i_state=1,N_states
! 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)

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

@ -14,7 +14,7 @@ BEGIN_PROVIDER [integer(bit_kind), psi_active, (N_int,2,psi_det_size)]
psi_active(j,1,i) = iand(psi_ref(j,1,i),cas_bitmask(j,1,1))
psi_active(j,2,i) = iand(psi_ref(j,2,i),cas_bitmask(j,1,1))
enddo
call debug_det(psi_active(1,1,i),N_int)
! call debug_det(psi_active(1,1,i),N_int)
enddo
END_PROVIDER
@ -330,6 +330,7 @@ subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final)
i_part = list_virt_reverse(p1)
do i_state = 1, N_states
delta_e_act(i_state) += one_creat_virt(i_hole,i_part,i_state)
! delta_e_act += 1.d12
enddo
else if (degree == 2)then
do i_state = 1, N_states

10
plugins/loc_cele/loc_cele.irp.f
View File

@ -101,10 +101,12 @@
cmoref = 0.d0
irot = 0
irot(1,1) = 5
irot(2,1) = 6
cmoref(6,1,1) = 1d0
cmoref(26,2,1) = 1d0
irot(1,1) = 14
irot(2,1) = 15
! cmoref(6,1,1) = 1.d0
! cmoref(26,2,1) = 1.d0
cmoref(36,1,1) = 1.d0
cmoref(56,2,1) = 1.d0
! !!! H2O
! irot(1,1) = 4

76
src/Determinants/H_apply.irp.f
View File

@ -195,6 +195,7 @@ subroutine copy_H_apply_buffer_to_wf
!call remove_duplicates_in_psi_det(found_duplicates)
end
subroutine remove_duplicates_in_psi_det(found_duplicates)
implicit none
logical, intent(out) :: found_duplicates
@ -270,6 +271,81 @@ subroutine remove_duplicates_in_psi_det(found_duplicates)
deallocate (duplicate,bit_tmp)
end
subroutine remove_duplicates_in_psi_det_new(found_duplicates)
implicit none
logical, intent(out) :: found_duplicates
BEGIN_DOC
! Removes duplicate determinants in the wave function.
END_DOC
integer :: i,j,k
integer(bit_kind), allocatable :: bit_tmp(:)
logical,allocatable :: duplicate(:)
allocate (duplicate(N_det), bit_tmp(N_det))
do i=1,N_det
integer, external :: det_search_key
!$DIR FORCEINLINE
bit_tmp(i) = det_search_key(psi_det_sorted_bit(1,1,i),N_int)
duplicate(i) = .False.
enddo
do i=1,N_det-1
if (duplicate(i)) then
cycle
endif
j = i+1
do while (bit_tmp(j)==bit_tmp(i))
if (duplicate(j)) then
j += 1
if (j > N_det) then
exit
else
cycle
endif
endif
duplicate(j) = .True.
do k=1,N_int
if ( (psi_det_sorted_bit(k,1,i) /= psi_det_sorted_bit(k,1,j) ) &
.or. (psi_det_sorted_bit(k,2,i) /= psi_det_sorted_bit(k,2,j) ) ) then
duplicate(j) = .False.
exit
endif
enddo
j += 1
if (j > N_det) then
exit
endif
enddo
enddo
found_duplicates = .False.
do i=1,N_det
if (duplicate(i)) then
found_duplicates = .True.
exit
endif
enddo
if (found_duplicates) then
k=0
do i=1,N_det
if (.not.duplicate(i)) then
k += 1
psi_det(:,:,k) = psi_det_sorted_bit (:,:,i)
psi_coef(k,:) = psi_coef_sorted_bit(i,:)
else
psi_det(:,:,k) = 0_bit_kind
psi_coef(k,:) = 0.d0
endif
enddo
N_det = k
call write_bool(output_determinants,found_duplicates,'Found duplicate determinants')
SOFT_TOUCH N_det psi_det psi_coef
endif
deallocate (duplicate,bit_tmp)
end
subroutine fill_H_apply_buffer_no_selection(n_selected,det_buffer,Nint,iproc)
use bitmasks

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

@ -134,7 +134,6 @@ END_PROVIDER
!$OMP END CRITICAL
deallocate(tmp_a,tmp_b)
!$OMP END PARALLEL
END_PROVIDER
BEGIN_PROVIDER [ double precision, one_body_single_double_dm_mo_alpha, (mo_tot_num_align,mo_tot_num) ]

46
src/Determinants/truncate_wf.irp.f
View File

@ -1,8 +1,52 @@
program s2_eig_restart
implicit none
read_wf = .True.
call routine
call routine_2
end
subroutine routine_2
implicit none
integer :: i,j,k,l
use bitmasks
integer :: n_det_restart,degree
integer(bit_kind),allocatable :: psi_det_tmp(:,:,:)
double precision ,allocatable :: psi_coef_tmp(:,:),accu(:)
integer, allocatable :: index_restart(:)
allocate(index_restart(N_det))
print*, 'How many Slater determinants would ou like ?'
read(5,*)N_det_restart
do i = 1, N_det_restart
index_restart(i) = i
enddo
allocate (psi_det_tmp(N_int,2,N_det_restart),psi_coef_tmp(N_det_restart,N_states),accu(N_states))
accu = 0.d0
do i = 1, N_det_restart
do j = 1, N_int
psi_det_tmp(j,1,i) = psi_det(j,1,index_restart(i))
psi_det_tmp(j,2,i) = psi_det(j,2,index_restart(i))
enddo
do j = 1,N_states
psi_coef_tmp(i,j) = psi_coef(index_restart(i),j)
accu(j) += psi_coef_tmp(i,j) * psi_coef_tmp(i,j)
enddo
enddo
do j = 1, N_states
accu(j) = 1.d0/dsqrt(accu(j))
enddo
do j = 1,N_states
do i = 1, N_det_restart
psi_coef_tmp(i,j) = psi_coef_tmp(i,j) * accu(j)
enddo
enddo
call save_wavefunction_general(N_det_restart,N_states,psi_det_tmp,N_det_restart,psi_coef_tmp)
deallocate (psi_det_tmp,psi_coef_tmp,accu,index_restart)
end
subroutine routine
implicit none
call make_s2_eigenfunction