LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-12-23 04:43:50 +01:00
Code Issues Releases Wiki Activity

Merged eginer-master

Browse Source
This commit is contained in:
Anthony Scemama 2015-03-19 21:14:52 +01:00
parent 39bf2889ce
commit 5c957cf1f2
80 changed files with 2714 additions and 274 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
data/ezfio_defaults
View File

@ -23,13 +23,20 @@ determinants
threshold_selectors 0.999
read_wf False
s2_eig False
only_single_double_dm False
full_ci
n_det_max_fci 10000
n_det_max_fci_property 50000
pt2_max 1.e-4
do_pt2_end True
var_pt2_ratio 0.75
all_singles
n_det_max_fci 50000
pt2_max 1.e-8
do_pt2_end False
cassd
n_det_max_cassd 10000
pt2_max 1.e-4
@ -39,8 +46,14 @@ hartree_fock
n_it_scf_max 200
thresh_scf 1.e-10
cisd_selected
n_det_max_cisd 10000
pt2_max 1.e-4
cisd_sc2_selected
n_det_max_cisd_sc2 10000
pt2_max 1.e-4
do_pt2_end True
properties
z_one_point 3.9

33
scripts/generate_h_apply.py
View File

@ -13,6 +13,8 @@ initialization
declarations
decls_main
keys_work
do_double_excitations
check_double_excitation
copy_buffer
finalization
generate_psi_guess
@ -23,7 +25,8 @@ deinit_thread
skip
init_main
filter_integrals
filter2h2p
filterhole
filterparticle
""".split()
class H_apply(object):
@ -116,12 +119,23 @@ class H_apply(object):
buffer = buffer.replace('$'+key, value)
return buffer
def set_filter_2h_2p(self):
self["filter2h2p"] = """
! ! DIR$ FORCEINLINE
if(is_a_two_holes_two_particles(key))cycle
def unset_double_excitations(self):
self["do_double_excitations"] = ".False."
self["check_double_excitation"] = """
check_double_excitation = .False.
"""
def set_filter_holes(self):
self["filterhole"] = """
if(iand(ibset(0_bit_kind,j),hole(k,other_spin)).eq.0_bit_kind ) then
cycle
endif
"""
def set_filter_particl(self):
self["filterparticle"] = """
if(iand(ibset(0_bit_kind,j_a),hole(k_a,other_spin)).eq.0_bit_kind ) then
cycle
endif
"""
def set_perturbation(self,pert):
if self.perturbation is not None:
@ -165,9 +179,14 @@ class H_apply(object):
PROVIDE CI_electronic_energy psi_selectors_coef psi_selectors E_corr_per_selectors psi_det_sorted_bit
"""
self.data["keys_work"] = """
if(check_double_excitation)then
call perturb_buffer_%s(i_generator,keys_out,key_idx,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert, &
sum_norm_pert,sum_H_pert_diag,N_st,N_int)
"""%(pert,)
else
call perturb_buffer_by_mono_%s(i_generator,keys_out,key_idx,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert, &
sum_norm_pert,sum_H_pert_diag,N_st,N_int)
endif
"""%(pert,pert)
self.data["finalization"] = """
"""
self.data["copy_buffer"] = ""

3
src/BiInts/ao_bi_integrals.irp.f
View File

@ -354,6 +354,7 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ]
real :: map_mb
if (read_ao_integrals) then
integer :: load_ao_integrals
print*,'Reading the AO integrals'
if (load_ao_integrals(trim(ezfio_filename)//'/work/ao_integrals.bin') == 0) then
write(output_BiInts,*) 'AO integrals provided'
ao_bielec_integrals_in_map = .True.
@ -373,7 +374,7 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ]
PROVIDE progress_bar
call omp_init_lock(lock)
lmax = ao_num*(ao_num+1)/2
write(output_BiInts,*) 'providing the AO integrals'
write(output_BiInts,*) 'Providing the AO integrals'
call wall_time(wall_0)
call wall_time(wall_1)
call cpu_time(cpu_1)

3
src/BiInts/map_integrals.irp.f
View File

@ -276,7 +276,6 @@ double precision function get_mo_bielec_integral(i,j,k,l,map)
implicit none
BEGIN_DOC
! Returns one integral <ij|kl> in the MO basis
! i(1)j(1) 1/r12 k(2)l(2)
END_DOC
integer, intent(in) :: i,j,k,l
integer*8 :: idx
@ -293,7 +292,6 @@ double precision function mo_bielec_integral(i,j,k,l)
implicit none
BEGIN_DOC
! Returns one integral <ij|kl> in the MO basis
! i(1)j(1) 1/r12 k(2)l(2)
END_DOC
integer, intent(in) :: i,j,k,l
double precision :: get_mo_bielec_integral
@ -308,7 +306,6 @@ subroutine get_mo_bielec_integrals(j,k,l,sze,out_val,map)
BEGIN_DOC
! Returns multiple integrals <ij|kl> in the MO basis, all
! i for j,k,l fixed.
! i(1)j(1) 1/r12 k(2)l(2)
END_DOC
integer, intent(in) :: j,k,l, sze
real(integral_kind), intent(out) :: out_val(sze)

1
src/BiInts/mo_bi_integrals.irp.f
View File

@ -28,6 +28,7 @@ BEGIN_PROVIDER [ logical, mo_bielec_integrals_in_map ]
mo_bielec_integrals_in_map = .True.
if (read_mo_integrals) then
integer :: load_mo_integrals
print*,'Reading the MO integrals'
if (load_mo_integrals(trim(ezfio_filename)//'/work/mo_integrals.bin') == 0) then
write(output_BiInts,*) 'MO integrals provided'
return

3
src/CAS_SD_selected/NEEDED_MODULES
View File

@ -1 +1,2 @@
AOs BiInts Bitmask Dets Electrons Ezfio_files Generators_CAS Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Selectors_full Utils
AOs BiInts Bitmask Dets Electrons Ezfio_files Generators_CAS Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Properties Selectors_full Utils

2
src/CID_SC2_selected/NEEDED_MODULES
View File

@ -1,2 +1,2 @@
AOs BiInts Bitmask CISD CISD_selected Dets Electrons Ezfio_files Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Selectors_full SingleRefMethod Utils
AOs BiInts Bitmask CISD CISD_selected Dets Electrons Ezfio_files Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Properties Selectors_full SingleRefMethod Utils

1
src/CID_SC2_selected/cid_sc2_selection.irp.f
View File

@ -68,6 +68,5 @@ program cisd_sc2_selected
print*,'degree of excitation of such determinant : ',degree
enddo
print*,'coucou'
deallocate(pt2,norm_pert,H_pert_diag)
end

2
src/CID_selected/NEEDED_MODULES
View File

@ -1,2 +1,2 @@
AOs BiInts Bitmask CISD Dets Electrons Ezfio_files Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Selectors_full SingleRefMethod Utils
AOs BiInts Bitmask CISD Dets Electrons Ezfio_files Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Properties Selectors_full SingleRefMethod Utils

4
src/CISD_SC2_selected/H_apply.irp.f
View File

@ -3,6 +3,10 @@ BEGIN_SHELL [ /usr/bin/env python ]
from generate_h_apply import *
from perturbation import perturbations
s = H_apply("SC2_selected",SingleRef=True)
s.set_selection_pt2("epstein_nesbet_sc2_no_projected")
print s
s = H_apply("PT2",SingleRef=True)
s.set_perturbation("epstein_nesbet_sc2_no_projected")
print s

2
src/CISD_SC2_selected/NEEDED_MODULES
View File

@ -1,2 +1,2 @@
AOs BiInts Bitmask CISD CISD_selected Dets Electrons Ezfio_files Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Selectors_full SingleRefMethod Utils
AOs BiInts Bitmask CISD CISD_selected Dets Electrons Ezfio_files Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Properties Selectors_full SingleRefMethod Utils

6
src/CISD_SC2_selected/cisd_sc2_selection.irp.f
View File

@ -37,8 +37,7 @@ program cisd_sc2_selected
do while (N_det < n_det_max_cisd_sc2.and.maxval(abs(pt2(1:N_st))) > pt2_max)
print*,'----'
print*,''
call H_apply_cisd_selection(perturbation,pt2, norm_pert, H_pert_diag, N_st)
! soft_touch det_connections
call H_apply_SC2_selected(pt2, norm_pert, H_pert_diag, N_st)
call diagonalize_CI_SC2
print *, 'N_det = ', N_det
do i = 1, N_st
@ -51,13 +50,14 @@ program cisd_sc2_selected
endif
E_old(i) = CI_SC2_energy(i)
enddo
! print *, 'E corr = ', (E_old(1)) - HF_energy
if(dabs(E_old(i) - CI_SC2_energy(i) ).le.1.d-12)then
i_count += 1
selection_criterion_factor = selection_criterion_factor * 0.5d0
if(i_count > 5)then
exit
endif
else
i_count = 0
endif
if (abort_all) then
exit

2
src/CISD_selected/NEEDED_MODULES
View File

@ -1,2 +1,2 @@
AOs BiInts Bitmask CISD Dets Electrons Ezfio_files Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Selectors_full SingleRefMethod Utils
AOs BiInts Bitmask CISD Dets Electrons Ezfio_files Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Properties Selectors_full SingleRefMethod Utils

3
src/CISD_selected/cisd.ezfio_config Normal file
View File

@ -0,0 +1,3 @@
cisd_selected
n_det_max_cisd integer
pt2_max double precision

8
src/CISD_selected/cisd_selection.irp.f
View File

@ -18,11 +18,13 @@ program cisd
print *, 'E = ', CI_energy(i)
enddo
E_old = CI_energy
do while (maxval(abs(pt2(1:N_st))) > 1.d-4)
do while (maxval(abs(pt2(1:N_st))) > pt2_max.and.n_det < n_det_max_cisd)
print*,'----'
print*,''
call H_apply_cisd_selection(perturbation,pt2, norm_pert, H_pert_diag, N_st)
call diagonalize_CI
psi_det = psi_det_sorted
psi_coef = psi_coef_sorted
print*,'N_det = ',N_det
do i = 1, N_st
print*,'state ',i
@ -36,5 +38,9 @@ program cisd
exit
endif
enddo
N_det = min(N_det,n_det_max_cisd)
touch N_det psi_det psi_coef
call diagonalize_CI
deallocate(pt2,norm_pert,H_pert_diag)
call save_wavefunction
end

34
src/CISD_selected/options.irp.f Normal file
View File

@ -0,0 +1,34 @@
BEGIN_PROVIDER [ integer, n_det_max_cisd ]
implicit none
BEGIN_DOC
! Get n_det_max_cisd from EZFIO file
END_DOC
logical :: has_n_det_max_cisd
PROVIDE ezfio_filename
call ezfio_has_cisd_selected_n_det_max_cisd(has_n_det_max_cisd)
if (has_n_det_max_cisd) then
call ezfio_get_cisd_selected_n_det_max_cisd(n_det_max_cisd)
else
n_det_max_cisd = 30000
call ezfio_set_cisd_selected_n_det_max_cisd(n_det_max_cisd)
endif
print*,'n_det_max_cisd = ',n_det_max_cisd
END_PROVIDER
BEGIN_PROVIDER [ double precision , pt2_max ]
implicit none
BEGIN_DOC
! Get pt2_max from EZFIO file
END_DOC
logical :: has_pt2_max
PROVIDE ezfio_filename
call ezfio_has_cisd_selected_pt2_max(has_pt2_max)
if (has_pt2_max) then
call ezfio_get_cisd_selected_pt2_max(pt2_max)
else
pt2_max = 1.d-9
call ezfio_set_cisd_selected_pt2_max(pt2_max)
endif
print*,'pt2_max = ',pt2_max
END_PROVIDER

3
src/DDCI_selected/NEEDED_MODULES
View File

@ -1 +1,2 @@
AOs BiInts Bitmask Dets Electrons Ezfio_files Generators_CAS Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Selectors_full Utils
AOs BiInts Bitmask Dets Electrons Ezfio_files Generators_CAS Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Properties Selectors_full Utils

18
src/Dets/H_apply_template.f
View File

@ -26,7 +26,6 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene
integer :: N_elec_in_key_hole_2(2),N_elec_in_key_part_2(2)
double precision :: mo_bielec_integral
logical :: is_a_two_holes_two_particles
integer, allocatable :: ia_ja_pairs(:,:,:)
integer, allocatable :: ib_jb_pairs(:,:)
double precision :: diag_H_mat_elem
@ -36,6 +35,11 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene
ifirst=1
endif
logical :: check_double_excitation
check_double_excitation = .True.
$initialization
$omp_parallel
@ -163,7 +167,6 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene
k = ishft(j_b-1,-bit_kind_shift)+1
l = j_b-ishft(k-1,bit_kind_shift)-1
key(k,other_spin) = ibset(key(k,other_spin),l)
$filter2h2p
key_idx += 1
do k=1,N_int
keys_out(k,1,key_idx) = key(k,1)
@ -212,7 +215,6 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene
k = ishft(j_b-1,-bit_kind_shift)+1
l = j_b-ishft(k-1,bit_kind_shift)-1
key(k,ispin) = ibset(key(k,ispin),l)
$filter2h2p
key_idx += 1
do k=1,N_int
keys_out(k,1,key_idx) = key(k,1)
@ -270,12 +272,17 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator,iproc $param
integer :: kk,pp,other_spin,key_idx
integer :: N_elec_in_key_hole_1(2),N_elec_in_key_part_1(2)
integer :: N_elec_in_key_hole_2(2),N_elec_in_key_part_2(2)
logical :: is_a_two_holes_two_particles
integer, allocatable :: ia_ja_pairs(:,:,:)
logical, allocatable :: array_pairs(:,:)
double precision :: diag_H_mat_elem
integer(omp_lock_kind), save :: lck, ifirst=0
logical :: check_double_excitation
check_double_excitation = .True.
$check_double_excitation
if (ifirst == 0) then
ifirst=1
!$ call omp_init_lock(lck)
@ -333,11 +340,12 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator,iproc $param
hole = key_in
k = ishft(i_a-1,-bit_kind_shift)+1
j = i_a-ishft(k-1,bit_kind_shift)-1
$filterhole
hole(k,ispin) = ibclr(hole(k,ispin),j)
k_a = ishft(j_a-1,-bit_kind_shift)+1
l_a = j_a-ishft(k_a-1,bit_kind_shift)-1
$filterparticle
hole(k_a,ispin) = ibset(hole(k_a,ispin),l_a)
$filter2h2p
key_idx += 1
do k=1,N_int
keys_out(k,1,key_idx) = hole(k,1)

108
src/Dets/connected_to_ref.irp.f
View File

@ -225,6 +225,114 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet)
end
integer function connected_to_ref_by_mono(key,keys,Nint,N_past_in,Ndet)
use bitmasks
implicit none
integer, intent(in) :: Nint, N_past_in, Ndet
integer(bit_kind), intent(in) :: keys(Nint,2,Ndet)
integer(bit_kind), intent(in) :: key(Nint,2)
integer :: N_past
integer :: i, l
integer :: degree_x2
logical :: det_is_not_or_may_be_in_ref, t
double precision :: hij_elec
! output : 0 : not connected
! i : connected to determinant i of the past
! -i : is the ith determinant of the refernce wf keys
ASSERT (Nint > 0)
ASSERT (Nint == N_int)
connected_to_ref_by_mono = 0
N_past = max(1,N_past_in)
if (Nint == 1) then
do i=N_past-1,1,-1
degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + &
popcnt(xor( key(1,2), keys(1,2,i)))
if (degree_x2 > 3.and. degree_x2 <5) then
cycle
else if (degree_x2 == 4)then
cycle
else if(degree_x2 == 2)then
connected_to_ref_by_mono = i
return
endif
enddo
return
else if (Nint==2) then
do i=N_past-1,1,-1
degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + &
popcnt(xor( key(1,2), keys(1,2,i))) + &
popcnt(xor( key(2,1), keys(2,1,i))) + &
popcnt(xor( key(2,2), keys(2,2,i)))
if (degree_x2 > 3.and. degree_x2 <5) then
cycle
else if (degree_x2 == 4)then
cycle
else if(degree_x2 == 2)then
connected_to_ref_by_mono = i
return
endif
enddo
return
else if (Nint==3) then
do i=N_past-1,1,-1
degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + &
popcnt(xor( key(1,2), keys(1,2,i))) + &
popcnt(xor( key(2,1), keys(2,1,i))) + &
popcnt(xor( key(2,2), keys(2,2,i))) + &
popcnt(xor( key(3,1), keys(3,1,i))) + &
popcnt(xor( key(3,2), keys(3,2,i)))
if (degree_x2 > 3.and. degree_x2 <5) then
cycle
else if (degree_x2 == 4)then
cycle
else if(degree_x2 == 2)then
connected_to_ref_by_mono = i
return
endif
enddo
return
else
do i=N_past-1,1,-1
degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + &
popcnt(xor( key(1,2), keys(1,2,i)))
!DEC$ LOOP COUNT MIN(3)
do l=2,Nint
degree_x2 = degree_x2 + popcnt(xor( key(l,1), keys(l,1,i))) +&
popcnt(xor( key(l,2), keys(l,2,i)))
enddo
if (degree_x2 > 3.and. degree_x2 <5) then
cycle
else if (degree_x2 == 4)then
cycle
else if(degree_x2 == 2)then
connected_to_ref_by_mono = i
return
endif
enddo
endif
end
logical function det_is_not_or_may_be_in_ref(key,Nint)
use bitmasks
implicit none

99
src/Dets/density_matrix.irp.f
View File

@ -13,6 +13,11 @@
integer :: exc(0:2,2,2),n_occ_alpha
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) &
@ -68,6 +73,92 @@
deallocate(tmp_a,tmp_b)
!$OMP BARRIER
!$OMP END PARALLEL
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, one_body_single_double_dm_mo_alpha, (mo_tot_num_align,mo_tot_num) ]
&BEGIN_PROVIDER [ double precision, one_body_single_double_dm_mo_beta, (mo_tot_num_align,mo_tot_num) ]
implicit none
BEGIN_DOC
! Alpha and beta one-body density matrix for each state
END_DOC
integer :: j,k,l,m
integer :: occ(N_int*bit_kind_size,2)
double precision :: ck, cl, ckl
double precision :: phase
integer :: h1,h2,p1,p2,s1,s2, degree
integer :: exc(0:2,2,2),n_occ_alpha
double precision, allocatable :: tmp_a(:,:), tmp_b(:,:)
integer :: degree_respect_to_HF_k
integer :: degree_respect_to_HF_l
PROVIDE elec_alpha_num elec_beta_num
one_body_single_double_dm_mo_alpha = 0.d0
one_body_single_double_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_alpha,degree_respect_to_HF_k,degree_respect_to_HF_l)&
!$OMP SHARED(ref_bitmask,psi_det,psi_coef,N_int,N_states,state_average_weight,elec_alpha_num,&
!$OMP elec_beta_num,one_body_single_double_dm_mo_alpha,one_body_single_double_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) )
tmp_a = 0.d0
tmp_b = 0.d0
!$OMP DO SCHEDULE(dynamic)
do k=1,N_det
call bitstring_to_list(psi_det(1,1,k), occ(1,1), n_occ_alpha, N_int)
call bitstring_to_list(psi_det(1,2,k), occ(1,2), n_occ_alpha, N_int)
call get_excitation_degree(ref_bitmask,psi_det(1,1,k),degree_respect_to_HF_k,N_int)
do m=1,N_states
ck = psi_coef(k,m)*psi_coef(k,m) * state_average_weight(m)
call get_excitation_degree(ref_bitmask,psi_det(1,1,k),degree_respect_to_HF_l,N_int)
if(degree_respect_to_HF_l.le.0)then
do l=1,elec_alpha_num
j = occ(l,1)
tmp_a(j,j) += ck
enddo
do l=1,elec_beta_num
j = occ(l,2)
tmp_b(j,j) += ck
enddo
endif
enddo
do l=1,k-1
call get_excitation_degree(ref_bitmask,psi_det(1,1,l),degree_respect_to_HF_l,N_int)
if(degree_respect_to_HF_k.ne.0)cycle
if(degree_respect_to_HF_l.eq.2.and.degree_respect_to_HF_k.ne.2)cycle
call get_excitation_degree(psi_det(1,1,k),psi_det(1,1,l),degree,N_int)
if (degree /= 1) then
cycle
endif
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)
if (s1==1) then
tmp_a(h1,p1) += ckl
tmp_a(p1,h1) += ckl
else
tmp_b(h1,p1) += ckl
tmp_b(p1,h1) += ckl
endif
enddo
enddo
enddo
!$OMP END DO NOWAIT
!$OMP CRITICAL
one_body_single_double_dm_mo_alpha = one_body_single_double_dm_mo_alpha + tmp_a
!$OMP END CRITICAL
!$OMP CRITICAL
one_body_single_double_dm_mo_beta = one_body_single_double_dm_mo_beta + tmp_b
!$OMP END CRITICAL
deallocate(tmp_a,tmp_b)
!$OMP BARRIER
!$OMP END PARALLEL
END_PROVIDER
BEGIN_PROVIDER [ double precision, one_body_dm_mo, (mo_tot_num_align,mo_tot_num) ]
@ -78,6 +169,14 @@ BEGIN_PROVIDER [ double precision, one_body_dm_mo, (mo_tot_num_align,mo_tot_num)
one_body_dm_mo = one_body_dm_mo_alpha + one_body_dm_mo_beta
END_PROVIDER
BEGIN_PROVIDER [ double precision, one_body_spin_density_mo, (mo_tot_num_align,mo_tot_num) ]
implicit none
BEGIN_DOC
! rho(alpha) - rho(beta)
END_DOC
one_body_spin_density_mo = one_body_dm_mo_alpha - one_body_dm_mo_beta
END_PROVIDER
subroutine set_natural_mos
implicit none
BEGIN_DOC

1
src/Dets/determinants.ezfio_config
View File

@ -16,4 +16,5 @@ determinants
read_wf logical
expected_s2 double precision
s2_eig logical
only_single_double_dm logical

36
src/Dets/determinants.irp.f
View File

@ -29,6 +29,20 @@ BEGIN_PROVIDER [ integer, N_det ]
ASSERT (N_det > 0)
END_PROVIDER
BEGIN_PROVIDER [integer, max_degree_exc]
implicit none
integer :: i,degree
max_degree_exc = 0
BEGIN_DOC
! Maximum degree of excitation in the wf
END_DOC
do i = 1, N_det
call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int)
if(degree.gt.max_degree_exc)then
max_degree_exc= degree
endif
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer, psi_det_size ]
implicit none
@ -729,6 +743,16 @@ subroutine save_wavefunction
call save_wavefunction_general(N_det,N_states,psi_det_sorted,size(psi_coef_sorted,1),psi_coef_sorted)
end
subroutine save_wavefunction_unsorted
implicit none
use bitmasks
BEGIN_DOC
! Save the wave function into the EZFIO file
END_DOC
call save_wavefunction_general(N_det,N_states,psi_det,size(psi_coef,1),psi_coef)
end
subroutine save_wavefunction_general(ndet,nstates,psidet,dim_psicoef,psicoef)
implicit none
BEGIN_DOC
@ -792,11 +816,23 @@ subroutine save_wavefunction_general(ndet,nstates,psidet,dim_psicoef,psicoef)
progress_bar(1) = 7
progress_value = dble(progress_bar(1))
allocate (psi_coef_save(ndet,nstates))
double precision :: accu_norm(nstates)
accu_norm = 0.d0
do k=1,nstates
do i=1,ndet
accu_norm(k) = accu_norm(k) + psicoef(i,k) * psicoef(i,k)
psi_coef_save(i,k) = psicoef(i,k)
enddo
enddo
do k = 1, nstates
accu_norm(k) = 1.d0/dsqrt(accu_norm(k))
enddo
do k=1,nstates
do i=1,ndet
psi_coef_save(i,k) = psi_coef_save(i,k) * accu_norm(k)
enddo
enddo
call ezfio_set_determinants_psi_coef(psi_coef_save)
call write_int(output_dets,ndet,'Saved determinants')
call stop_progress

2
src/Dets/diagonalize_CI.irp.f
View File

@ -69,10 +69,8 @@ END_PROVIDER
i_state = 0
do j=1,N_det
call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2)
! print *, 'j = ',j,s2, expected_s2
if(dabs(s2-expected_s2).le.0.3d0)then
i_state += 1
! print *, 'i_state = ',i_state
do i=1,N_det
CI_eigenvectors(i,i_state) = eigenvectors(i,j)
enddo

3
src/Dets/diagonalize_CI_SC2.irp.f
View File

@ -35,7 +35,8 @@ END_PROVIDER
do i=1,N_det
CI_SC2_eigenvectors(i,j) = psi_coef(i,j)
enddo
CI_SC2_electronic_energy(j) = CI_electronic_energy(j)
! TODO : check comment
! CI_SC2_electronic_energy(j) = CI_electronic_energy(j)
enddo
call CISD_SC2(psi_det,CI_SC2_eigenvectors,CI_SC2_electronic_energy, &

72
src/Dets/diagonalize_CI_mono.irp.f Normal file
View File

@ -0,0 +1,72 @@
BEGIN_PROVIDER [ double precision, CI_electronic_energy_mono, (N_states_diag) ]
&BEGIN_PROVIDER [ double precision, CI_eigenvectors_mono, (N_det,N_states_diag) ]
&BEGIN_PROVIDER [ double precision, CI_eigenvectors_s2_mono, (N_states_diag) ]
implicit none
BEGIN_DOC
! Eigenvectors/values of the CI matrix
END_DOC
integer :: i,j
do j=1,N_states_diag
do i=1,N_det
CI_eigenvectors_mono(i,j) = psi_coef(i,j)
enddo
enddo
if (diag_algorithm == "Davidson") then
call davidson_diag(psi_det,CI_eigenvectors_mono,CI_electronic_energy, &
size(CI_eigenvectors_mono,1),N_det,N_states_diag,N_int,output_Dets)
else if (diag_algorithm == "Lapack") then
double precision, allocatable :: eigenvectors(:,:), eigenvalues(:)
allocate (eigenvectors(size(H_matrix_all_dets,1),N_det))
allocate (eigenvalues(N_det))
call lapack_diag(eigenvalues,eigenvectors, &
H_matrix_all_dets,size(H_matrix_all_dets,1),N_det)
CI_electronic_energy_mono(:) = 0.d0
do i=1,N_det
CI_eigenvectors_mono(i,1) = eigenvectors(i,1)
enddo
integer :: i_state
double precision :: s2
i_state = 0
do j=1,N_det
call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2)
if(dabs(s2-expected_s2).le.0.3d0)then
print*,'j = ',j
print*,'e = ',eigenvalues(j)
print*,'c = ',dabs(eigenvectors(1,j))
if(dabs(eigenvectors(1,j)).gt.0.9d0)then
i_state += 1
do i=1,N_det
CI_eigenvectors_mono(i,i_state) = eigenvectors(i,j)
enddo
CI_electronic_energy_mono(i_state) = eigenvalues(j)
CI_eigenvectors_s2_mono(i_state) = s2
endif
endif
if (i_state.ge.N_states_diag) then
exit
endif
enddo
deallocate(eigenvectors,eigenvalues)
endif
END_PROVIDER
subroutine diagonalize_CI_mono
implicit none
BEGIN_DOC
! Replace the coefficients of the CI states by the coefficients of the
! eigenstates of the CI matrix
END_DOC
integer :: i,j
do j=1,N_states_diag
do i=1,N_det
psi_coef(i,j) = CI_eigenvectors_mono(i,j)
enddo
enddo
SOFT_TOUCH psi_coef CI_electronic_energy_mono CI_eigenvectors_mono CI_eigenvectors_s2_mono
end

16
src/Dets/excitations_utils.irp.f Normal file
View File

@ -0,0 +1,16 @@
subroutine apply_mono(i_hole,i_particle,ispin_excit,key_in,Nint)
implicit none
integer, intent(in) :: i_hole,i_particle,ispin_excit,Nint
integer(bit_kind), intent(inout) :: key_in(Nint,2)
integer :: k,j
use bitmasks
! hole
k = ishft(i_hole-1,-bit_kind_shift)+1
j = i_hole-ishft(k-1,bit_kind_shift)-1
key_in(k,ispin_excit) = ibclr(key_in(k,ispin_excit),j)
k = ishft(i_particle-1,-bit_kind_shift)+1
j = i_particle-ishft(k-1,bit_kind_shift)-1
key_in(k,ispin_excit) = ibset(key_in(k,ispin_excit),j)
end

79
src/Dets/guess_doublet.irp.f Normal file
View File

@ -0,0 +1,79 @@
program put_gess
use bitmasks
implicit none
integer :: i,j,N_det_tmp,N_states_tmp
integer :: list(N_int*bit_kind_size,2)
integer(bit_kind) :: string(N_int,2)
integer(bit_kind) :: psi_det_tmp(N_int,2,3)
double precision :: psi_coef_tmp(3,1)
integer :: iorb,jorb,korb
print*,'which open shells ?'
read(5,*)iorb,jorb,korb
print*,iorb,jorb,korb
N_states= 1
N_det= 3
list = 0
list(1,1) = 1
list(1,2) = 1
list(2,1) = 2
list(2,2) = 2
list(3,1) = iorb
list(4,1) = jorb
list(3,2) = korb
print*,'passed'
call list_to_bitstring( string(1,1), list(1,1), elec_alpha_num, N_int)
print*,'passed'
call list_to_bitstring( string(1,2), list(1,2), elec_beta_num, N_int)
print*,'passed'
call print_det(string,N_int)
do j = 1,2
do i = 1, N_int
psi_det(i,j,1) = string(i,j)
enddo
enddo
psi_coef(1,1) = 1.d0/dsqrt(3.d0)
print*,'passed 1'
list = 0
list(1,1) = 1
list(1,2) = 1
list(2,1) = 2
list(2,2) = 2
list(3,1) = iorb
list(4,1) = korb
list(3,2) = jorb
call list_to_bitstring( string(1,1), list(1,1), elec_alpha_num, N_int)
call list_to_bitstring( string(1,2), list(1,2), elec_beta_num, N_int)
call print_det(string,N_int)
do j = 1,2
do i = 1, N_int
psi_det(i,j,2) = string(i,j)
enddo
enddo
psi_coef(2,1) = 1.d0/dsqrt(3.d0)
print*,'passed 2'
list = 0
list(1,1) = 1
list(1,2) = 1
list(2,1) = 2
list(2,2) = 2
list(3,1) = korb
list(4,1) = jorb
list(3,2) = iorb
call list_to_bitstring( string(1,1), list(1,1), elec_alpha_num, N_int)
call list_to_bitstring( string(1,2), list(1,2), elec_beta_num, N_int)
call print_det(string,N_int)
do j = 1,2
do i = 1, N_int
psi_det(i,j,3) = string(i,j)
enddo
enddo
psi_coef(3,1) = 1.d0/dsqrt(3.d0)
print*,'passed 3'
call save_wavefunction
end

44
src/Dets/guess_singlet.irp.f Normal file
View File

@ -0,0 +1,44 @@
program put_gess
use bitmasks
implicit none
integer :: i,j,N_det_tmp,N_states_tmp
integer :: list(N_int*bit_kind_size,2)
integer(bit_kind) :: string(N_int,2)
integer(bit_kind) :: psi_det_tmp(N_int,2,2)
double precision :: psi_coef_tmp(2,1)
integer :: iorb,jorb
print*,'which open shells ?'
read(5,*)iorb,jorb
N_states= 1
N_det= 2
list = 0
list(1,1) = iorb
list(1,2) = jorb
call list_to_bitstring( string(1,1), list(1,1), elec_alpha_num, N_int)
call list_to_bitstring( string(1,2), list(1,2), elec_beta_num, N_int)
call print_det(string,N_int)
do j = 1,2
do i = 1, N_int
psi_det(i,j,1) = string(i,j)
enddo
enddo
psi_coef(1,1) = 1.d0/dsqrt(2.d0)
list = 0
list(1,1) = jorb
list(1,2) = iorb
call list_to_bitstring( string(1,1), list(1,1), elec_alpha_num, N_int)
call list_to_bitstring( string(1,2), list(1,2), elec_beta_num, N_int)
call print_det(string,N_int)
do j = 1,2
do i = 1, N_int
psi_det(i,j,2) = string(i,j)
enddo
enddo
psi_coef(2,1) = 1.d0/dsqrt(2.d0)
call save_wavefunction
end

48
src/Dets/guess_triplet.irp.f Normal file
View File

@ -0,0 +1,48 @@
program put_gess
use bitmasks
implicit none
integer :: i,j,N_det_tmp,N_states_tmp
integer :: list(N_int*bit_kind_size,2)
integer(bit_kind) :: string(N_int,2)
integer(bit_kind) :: psi_det_tmp(N_int,2,2)
double precision :: psi_coef_tmp(2,1)
integer :: iorb,jorb
print*,'which open shells ?'
read(5,*)iorb,jorb
N_states= 1
N_det= 2
print*,'iorb = ',iorb
print*,'jorb = ',jorb
list = 0
list(1,1) = iorb
list(1,2) = jorb
string = 0
call list_to_bitstring( string(1,1), list(1,1), elec_alpha_num, N_int)
call list_to_bitstring( string(1,2), list(1,2), elec_beta_num, N_int)
call print_det(string,N_int)
do j = 1,2
do i = 1, N_int
psi_det(i,j,1) = string(i,j)
enddo
enddo
psi_coef(1,1) = 1.d0/dsqrt(2.d0)
list = 0
list(1,1) = jorb
list(1,2) = iorb
string = 0
call list_to_bitstring( string(1,1), list(1,1), elec_alpha_num, N_int)
call list_to_bitstring( string(1,2), list(1,2), elec_beta_num, N_int)
call print_det(string,N_int)
do j = 1,2
do i = 1, N_int
psi_det(i,j,2) = string(i,j)
enddo
enddo
psi_coef(2,1) = -1.d0/dsqrt(2.d0)
call save_wavefunction
end

8
src/Dets/options.irp.f
View File

@ -22,6 +22,14 @@ T.set_ezfio_name( "read_wf" )
T.set_output ( "output_dets" )
print T
T.set_type ( "logical" )
T.set_name ( "only_single_double_dm" )
T.set_doc ( "If true, The One body DM is calculated with ignoring the Double<->Doubles extra diag elements" )
T.set_ezfio_name( "only_single_double_dm" )
T.set_output ( "output_dets" )
print T
T.set_name ( "s2_eig" )
T.set_doc ( "Force the wave function to be an eigenfunction of S^2" )
T.set_ezfio_name( "s2_eig" )

186
src/Dets/slater_rules.irp.f
View File

@ -488,6 +488,146 @@ end
subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble)
use bitmasks
implicit none
BEGIN_DOC
! Returns <i|H|j> where i and j are determinants
END_DOC
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2)
double precision, intent(out) :: hij,hmono,hdouble
integer :: exc(0:2,2,2)
integer :: degree
double precision :: get_mo_bielec_integral
integer :: m,n,p,q
integer :: i,j,k
integer :: occ(Nint*bit_kind_size,2)
double precision :: diag_H_mat_elem, phase,phase_2
integer :: n_occ_alpha, n_occ_beta
logical :: has_mipi(Nint*bit_kind_size)
double precision :: mipi(Nint*bit_kind_size), miip(Nint*bit_kind_size)
PROVIDE mo_bielec_integrals_in_map mo_integrals_map
ASSERT (Nint > 0)
ASSERT (Nint == N_int)
ASSERT (sum(popcnt(key_i(:,1))) == elec_alpha_num)
ASSERT (sum(popcnt(key_i(:,2))) == elec_beta_num)
ASSERT (sum(popcnt(key_j(:,1))) == elec_alpha_num)
ASSERT (sum(popcnt(key_j(:,2))) == elec_beta_num)
hij = 0.d0
hmono = 0.d0
hdouble = 0.d0
!DEC$ FORCEINLINE
call get_excitation_degree(key_i,key_j,degree,Nint)
select case (degree)
case (2)
call get_double_excitation(key_i,key_j,exc,phase,Nint)
if (exc(0,1,1) == 1) then
! Mono alpha, mono beta
hij = phase*get_mo_bielec_integral( &
exc(1,1,1), &
exc(1,1,2), &
exc(1,2,1), &
exc(1,2,2) ,mo_integrals_map)
else if (exc(0,1,1) == 2) then
! Double alpha
hij = phase*(get_mo_bielec_integral( &
exc(1,1,1), &
exc(2,1,1), &
exc(1,2,1), &
exc(2,2,1) ,mo_integrals_map) - &
get_mo_bielec_integral( &
exc(1,1,1), &
exc(2,1,1), &
exc(2,2,1), &
exc(1,2,1) ,mo_integrals_map) )
else if (exc(0,1,2) == 2) then
! Double beta
hij = phase*(get_mo_bielec_integral( &
exc(1,1,2), &
exc(2,1,2), &
exc(1,2,2), &
exc(2,2,2) ,mo_integrals_map) - &
get_mo_bielec_integral( &
exc(1,1,2), &
exc(2,1,2), &
exc(2,2,2), &
exc(1,2,2) ,mo_integrals_map) )
endif
case (1)
call get_mono_excitation(key_i,key_j,exc,phase,Nint)
call bitstring_to_list(key_i(1,1), occ(1,1), n_occ_alpha, Nint)
call bitstring_to_list(key_i(1,2), occ(1,2), n_occ_beta, Nint)
has_mipi = .False.
if (exc(0,1,1) == 1) then
! Mono alpha
m = exc(1,1,1)
p = exc(1,2,1)
do k = 1, elec_alpha_num
i = occ(k,1)
if (.not.has_mipi(i)) then
mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
has_mipi(i) = .True.
endif
enddo
do k = 1, elec_beta_num
i = occ(k,2)
if (.not.has_mipi(i)) then
mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
has_mipi(i) = .True.
endif
enddo
do k = 1, elec_alpha_num
hdouble = hdouble + mipi(occ(k,1)) - miip(occ(k,1))
enddo
do k = 1, elec_beta_num
hdouble = hdouble + mipi(occ(k,2))
enddo
else
! Mono beta
m = exc(1,1,2)
p = exc(1,2,2)
do k = 1, elec_beta_num
i = occ(k,2)
if (.not.has_mipi(i)) then
mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
has_mipi(i) = .True.
endif
enddo
do k = 1, elec_alpha_num
i = occ(k,1)
if (.not.has_mipi(i)) then
mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
has_mipi(i) = .True.
endif
enddo
do k = 1, elec_alpha_num
hdouble = hdouble + mipi(occ(k,1))
enddo
do k = 1, elec_beta_num
hdouble = hdouble + mipi(occ(k,2)) - miip(occ(k,2))
enddo
endif
hmono = mo_mono_elec_integral(m,p)
hij = phase*(hdouble + hmono)
case (0)
hij = diag_H_mat_elem(key_i,Nint)
end select
end
subroutine i_H_psi(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array)
use bitmasks
implicit none
@ -523,6 +663,52 @@ subroutine i_H_psi(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array)
enddo
end
subroutine i_H_psi_sec_ord(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array,idx_interaction,interactions)
use bitmasks
implicit none
integer, intent(in) :: Nint, Ndet,Ndet_max,Nstate
integer(bit_kind), intent(in) :: keys(Nint,2,Ndet)
integer(bit_kind), intent(in) :: key(Nint,2)
double precision, intent(in) :: coef(Ndet_max,Nstate)
double precision, intent(out) :: i_H_psi_array(Nstate)
double precision, intent(out) :: interactions(Ndet)
integer,intent(out) :: idx_interaction(0:Ndet)
integer :: i, ii,j
double precision :: phase
integer :: exc(0:2,2,2)
double precision :: hij
integer :: idx(0:Ndet),n_interact
BEGIN_DOC
! <key|H|psi> for the various Nstates
END_DOC
ASSERT (Nint > 0)
ASSERT (N_int == Nint)
ASSERT (Nstate > 0)
ASSERT (Ndet > 0)
ASSERT (Ndet_max >= Ndet)
i_H_psi_array = 0.d0
call filter_connected_i_H_psi0(keys,key,Nint,Ndet,idx)
n_interact = 0
do ii=1,idx(0)
i = idx(ii)
!DEC$ FORCEINLINE
call i_H_j(keys(1,1,i),key,Nint,hij)
if(dabs(hij).ge.1.d-8)then
if(i.ne.1)then
n_interact += 1
interactions(n_interact) = hij
idx_interaction(n_interact) = i
endif
endif
do j = 1, Nstate
i_H_psi_array(j) = i_H_psi_array(j) + coef(i,j)*hij
enddo
enddo
idx_interaction(0) = n_interact
end
subroutine i_H_psi_SC2(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array,idx_repeat)
use bitmasks

27
src/Full_CI/H_apply.irp.f
View File

@ -10,5 +10,32 @@ s = H_apply("FCI_PT2")
s.set_perturbation("epstein_nesbet_2x2")
print s
s = H_apply("FCI_mono")
s.set_selection_pt2("epstein_nesbet_2x2")
s.unset_double_excitations()
print s
s = H_apply("select_mono_delta_rho")
s.unset_double_excitations()
s.set_selection_pt2("delta_rho_one_point")
print s
s = H_apply("select_mono_di_delta_rho")
s.set_selection_pt2("delta_rho_one_point")
print s
s = H_apply("pt2_mono_delta_rho")
s.unset_double_excitations()
s.set_perturbation("delta_rho_one_point")
print s
s = H_apply("pt2_mono_di_delta_rho")
s.set_perturbation("delta_rho_one_point")
print s
END_SHELL

2
src/Full_CI/NEEDED_MODULES
View File

@ -1,2 +1,2 @@
AOs BiInts Bitmask Dets Electrons Ezfio_files Generators_full Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Selectors_full Utils
AOs BiInts Bitmask Dets Electrons Ezfio_files Generators_full Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Properties Selectors_full Utils

1
src/Full_CI/full_ci.ezfio_config
View File

@ -1,5 +1,6 @@
full_ci
n_det_max_fci integer
n_det_max_fci_property integer
pt2_max double precision
do_pt2_end logical
energy double precision

24
src/Full_CI/full_ci.irp.f
View File

@ -27,8 +27,19 @@ program full_ci
print *, 'E+PT2 = ', CI_energy+pt2
print *, '-----'
endif
double precision :: i_H_psi_array(N_states),diag_H_mat_elem,h,i_O1_psi_array(N_states)
if(read_wf)then
call i_H_psi(psi_det(1,1,N_det),psi_det,psi_coef,N_int,N_det,psi_det_size,N_states,i_H_psi_array)
h = diag_H_mat_elem(psi_det(1,1,N_det),N_int)
selection_criterion = dabs(psi_coef(N_det,1) * (i_H_psi_array(1) - h * psi_coef(N_det,1))) * 0.1d0
soft_touch selection_criterion
endif
integer :: n_det_before
print*,'Beginning the selection ...'
do while (N_det < n_det_max_fci.and.maxval(abs(pt2(1:N_st))) > pt2_max)
n_det_before = N_det
call H_apply_FCI(pt2, norm_pert, H_pert_diag, N_st)
PROVIDE psi_coef
@ -43,6 +54,9 @@ program full_ci
endif
call diagonalize_CI
call save_wavefunction
if(n_det_before == N_det)then
selection_criterion = selection_criterion * 0.5d0
endif
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
print *, 'PT2 = ', pt2
@ -55,16 +69,15 @@ program full_ci
endif
enddo
N_det = min(n_det_max_fci,N_det)
if(do_pt2_end)then
threshold_selectors = 1.d0
threshold_generators = 0.999d0
touch N_det psi_det psi_coef
call diagonalize_CI
if(do_pt2_end)then
print*,'Last iteration only to compute the PT2'
threshold_selectors = 1.d0
threshold_generators = 0.999d0
call H_apply_FCI_PT2(pt2, norm_pert, H_pert_diag, N_st)
print *, 'Final step'
!! call remove_small_contributions
!! call diagonalize_CI
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
print *, 'PT2 = ', pt2
@ -73,5 +86,6 @@ program full_ci
print *, '-----'
call ezfio_set_full_ci_energy_pt2(CI_energy+pt2)
endif
call save_wavefunction
deallocate(pt2,norm_pert)
end

8
src/Full_CI/options.irp.f
View File

@ -9,6 +9,14 @@ T.set_ezfio_name( "N_det_max_fci" )
T.set_output ( "output_full_ci" )
print T
T.set_type ( "integer" )
T.set_name ( "N_det_max_fci_property" )
T.set_doc ( "Max number of determinants in the wave function when you select for a given property" )
T.set_ezfio_dir ( "full_ci" )
T.set_ezfio_name( "N_det_max_fci_property" )
T.set_output ( "output_full_ci" )
print T
T.set_type ( "logical" )
T.set_name ( "do_pt2_end" )
T.set_doc ( "If true, compute the PT2 at the end of the selection" )

15
src/Generators_full/generators.irp.f
View File

@ -50,6 +50,21 @@ BEGIN_PROVIDER [ integer(bit_kind), psi_generators, (N_int,2,psi_det_size) ]
END_PROVIDER
BEGIN_PROVIDER [integer, degree_max_generators]
implicit none
BEGIN_DOC
! Max degree of excitation (respect to HF) of the generators
END_DOC
integer :: i,degree
degree_max_generators = 0
do i = 1, N_det_generators
call get_excitation_degree(HF_bitmask,psi_generators(1,1,i),degree,N_int)
if(degree .gt. degree_max_generators)then
degree_max_generators = degree
endif
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer, size_select_max]
implicit none
BEGIN_DOC

0
src/Primitive_basis/ASSUMPTIONS.rst → src/Generators_restart/ASSUMPTIONS.rst
View File

8
src/Generators_restart/Makefile Normal file
View File

@ -0,0 +1,8 @@
default: all
# Define here all new external source files and objects.Don't forget to prefix the
# object files with IRPF90_temp/
SRC=
OBJ=
include $(QPACKAGE_ROOT)/src/Makefile.common

1
src/Generators_restart/NEEDED_MODULES Normal file
View File

@ -0,0 +1 @@
AOs BiInts Bitmask Dets Electrons Ezfio_files MonoInts MOs Nuclei Output Utils

4
src/Generators_restart/README.rst Normal file
View File

@ -0,0 +1,4 @@
=========================
Generators_restart Module
=========================

59
src/Generators_restart/generators.irp.f Normal file
View File

@ -0,0 +1,59 @@
use bitmasks
BEGIN_PROVIDER [ integer, N_det_generators ]
implicit none
BEGIN_DOC
! Read the wave function
END_DOC
integer :: i
integer, save :: ifirst = 0
double precision :: norm
read_wf = .True.
if(ifirst == 0)then
N_det_generators = N_det
ifirst = 1
endif
call write_int(output_dets,N_det_generators,'Number of generators')
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), psi_generators, (N_int,2,N_det_generators) ]
&BEGIN_PROVIDER [ double precision, psi_generators_coef, (N_det_generators,N_states)
implicit none
BEGIN_DOC
! read wf
!
END_DOC
integer :: i, k
integer, save :: ifirst = 0
if(ifirst == 0)then
do i=1,N_det_generators
do k=1,N_int
psi_generators(k,1,i) = psi_det(k,1,i)
psi_generators(k,2,i) = psi_det(k,2,i)
enddo
do k = 1, N_states
psi_generators_coef(i,k) = psi_coef(i,k)
enddo
enddo
ifirst = 1
endif
END_PROVIDER
BEGIN_PROVIDER [ integer, size_select_max]
implicit none
BEGIN_DOC
! Size of the select_max array
END_DOC
size_select_max = 10000
END_PROVIDER
BEGIN_PROVIDER [ double precision, select_max, (size_select_max) ]
implicit none
BEGIN_DOC
! Memo to skip useless selectors
END_DOC
select_max = huge(1.d0)
END_PROVIDER

41
src/MOs/utils.irp.f
View File

@ -139,3 +139,44 @@ subroutine mo_as_eigvectors_of_mo_matrix_sort_by_observable(matrix,observable,n,
mo_label = label
SOFT_TOUCH mo_coef mo_label
end
subroutine mo_sort_by_observable(observable,label)
implicit none
character*(64), intent(in) :: label
double precision, intent(in) :: observable(mo_tot_num)
double precision, allocatable :: mo_coef_new(:,:),value(:)
integer,allocatable :: iorder(:)
allocate(mo_coef_new(ao_num_align,mo_tot_num),value(mo_tot_num),iorder(mo_tot_num))
print*,'allocate !'
mo_coef_new = mo_coef
do i = 1, mo_tot_num
iorder(i) = i
value(i) = observable(i)
enddo
integer :: i,j,k,index
print*,'sort ....'
call dsort(value,iorder,mo_tot_num)
do i = 1, mo_tot_num
index = iorder(i)
do j = 1, mo_tot_num
mo_coef(j,i) = mo_coef_new(j,index)
enddo
enddo
write (output_mos,'(A)'), 'MOs are now **'//trim(label)//'**'
write (output_mos,'(A)'), ''
deallocate(mo_coef_new,value)
! call write_time(output_mos)
mo_label = label
SOFT_TOUCH mo_coef mo_label
end

2
src/MP2/NEEDED_MODULES
View File

@ -1,2 +1,2 @@
AOs BiInts Bitmask Dets Electrons Ezfio_files Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Selectors_full SingleRefMethod Utils
AOs BiInts Bitmask Dets Electrons Ezfio_files Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Properties Selectors_full SingleRefMethod Utils

4
src/Molden/print_mo.irp.f
View File

@ -91,7 +91,11 @@ subroutine write_Ao_basis(i_unit_output)
! write(i_unit_output,*),j,i_shell,i_ao!trim(character_shell)
do k = 1, ao_prim_num(i_ao)
i_prim +=1
if(i_prim.lt.100)then
write(i_unit_output,'(4X,I3,3X,A1,6X,I2,6X,F16.7,2X,F16.12)')i_shell,character_shell,i_prim,ao_expo_unsorted(i_ao,k),ao_coef_unnormalized(i_ao,k)
else
write(i_unit_output,'(4X,I3,3X,A1,5X,I3,6X,F16.7,2X,F16.12)')i_shell,character_shell,i_prim,ao_expo_unsorted(i_ao,k),ao_coef_unnormalized(i_ao,k)
endif
enddo
write(i_unit_output,*)''
enddo

68
src/MonoInts/pot_ao_ints.irp.f
View File

@ -64,6 +64,73 @@
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral_per_atom, (ao_num_align,ao_num,nucl_num)]
BEGIN_DOC
! ao_nucl_elec_integral_per_atom(i,j,k) = -<AO(i)|1/|r-Rk|AO(j)>
! where Rk is the geometry of the kth atom
END_DOC
implicit none
double precision :: alpha, beta, gama, delta
integer :: i_c,num_A,num_B
double precision :: A_center(3),B_center(3),C_center(3)
integer :: power_A(3),power_B(3)
integer :: i,j,k,l,n_pt_in,m
double precision ::overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult
integer :: nucl_numC
! Important for OpenMP
ao_nucl_elec_integral_per_atom = 0.d0
do k = 1, nucl_num
C_center(1) = nucl_coord(k,1)
C_center(2) = nucl_coord(k,2)
C_center(3) = nucl_coord(k,3)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,l,m,alpha,beta,A_center,B_center,power_A,power_B, &
!$OMP num_A,num_B,c,n_pt_in) &
!$OMP SHARED (k,ao_num,ao_prim_num,ao_expo_transp,ao_power,ao_nucl,nucl_coord,ao_coef_transp, &
!$OMP n_pt_max_integrals,ao_nucl_elec_integral_per_atom,nucl_num,C_center)
n_pt_in = n_pt_max_integrals
!$OMP DO SCHEDULE (guided)
double precision :: c
do j = 1, ao_num
power_A(1)= ao_power(j,1)
power_A(2)= ao_power(j,2)
power_A(3)= ao_power(j,3)
num_A = ao_nucl(j)
A_center(1) = nucl_coord(num_A,1)
A_center(2) = nucl_coord(num_A,2)
A_center(3) = nucl_coord(num_A,3)
do i = 1, ao_num
power_B(1)= ao_power(i,1)
power_B(2)= ao_power(i,2)
power_B(3)= ao_power(i,3)
num_B = ao_nucl(i)
B_center(1) = nucl_coord(num_B,1)
B_center(2) = nucl_coord(num_B,2)
B_center(3) = nucl_coord(num_B,3)
c = 0.d0
do l=1,ao_prim_num(j)
alpha = ao_expo_transp(l,j)
do m=1,ao_prim_num(i)
beta = ao_expo_transp(m,i)
c = c + NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in) &
* ao_coef_transp(l,j)*ao_coef_transp(m,i)
enddo
enddo
ao_nucl_elec_integral_per_atom(i,j,k) = -c
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
enddo
END_PROVIDER
double precision function NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in)
! function that calculate the folowing integral :
@ -290,7 +357,6 @@ recursive subroutine I_x1_pol_mult_mono_elec(a,c,R1x,R1xp,R2x,d,nd,n_pt_in)
! print*,'nd_in = ',nd
if( (a==0) .and. (c==0))then
! print*,'coucou !'
nd = 0
d(0) = 1.d0
return

39
src/MonoInts/pot_mo_ints.irp.f
View File

@ -2,6 +2,9 @@ BEGIN_PROVIDER [double precision, mo_nucl_elec_integral, (mo_tot_num_align,mo_to
implicit none
integer :: i1,j1,i,j
double precision :: c_i1,c_j1
BEGIN_DOC
! interaction nuclear electron on the MO basis
END_DOC
mo_nucl_elec_integral = 0.d0
!$OMP PARALLEL DO DEFAULT(none) &
@ -11,9 +14,9 @@ BEGIN_PROVIDER [double precision, mo_nucl_elec_integral, (mo_tot_num_align,mo_to
do i = 1, mo_tot_num
do j = 1, mo_tot_num
do i1 = 1,ao_num
c_i1 = mo_coef(i1,i) ! <AO(i1)|MO(i)>
c_i1 = mo_coef(i1,i)
do j1 = 1,ao_num
c_j1 = c_i1*mo_coef(j1,j) ! <AO(j1)|MO(j)>
c_j1 = c_i1*mo_coef(j1,j)
mo_nucl_elec_integral(j,i) = mo_nucl_elec_integral(j,i) + &
c_j1 * ao_nucl_elec_integral(j1,i1)
enddo
@ -23,3 +26,35 @@ BEGIN_PROVIDER [double precision, mo_nucl_elec_integral, (mo_tot_num_align,mo_to
!$OMP END PARALLEL DO
END_PROVIDER
BEGIN_PROVIDER [double precision, mo_nucl_elec_integral_per_atom, (mo_tot_num_align,mo_tot_num,nucl_num)]
implicit none
integer :: i1,j1,i,j,k
double precision :: c_i1,c_j1
BEGIN_DOC
! mo_nucl_elec_integral_per_atom(i,j,k) = -<MO(i)|1/|r-Rk|MO(j)>
! where Rk is the geometry of the kth atom
END_DOC
mo_nucl_elec_integral_per_atom = 0.d0
do k = 1, nucl_num
!$OMP PARALLEL DO DEFAULT(none) &
!$OMP PRIVATE(i,j,i1,j1,c_j1,c_i1) &
!$OMP SHARED(mo_tot_num,ao_num,mo_coef, &
!$OMP mo_nucl_elec_integral_per_atom, ao_nucl_elec_integral_per_atom,k)
do i = 1, mo_tot_num
do j = 1, mo_tot_num
do i1 = 1,ao_num
c_i1 = mo_coef(i1,i)
do j1 = 1,ao_num
c_j1 = c_i1*mo_coef(j1,j)
mo_nucl_elec_integral_per_atom(j,i,k) = mo_nucl_elec_integral_per_atom(j,i,k) + &
c_j1 * ao_nucl_elec_integral_per_atom(j1,i1,k)
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
enddo
END_PROVIDER

6
src/MonoInts/spread_dipole_ao.irp.f
View File

@ -339,6 +339,12 @@ end
dist = (A_center - B_center)*(A_center - B_center)
P_center = (alpha * A_center + beta * B_center) * p_inv
factor = dexp(-rho * dist)
if(power_B == 0 .and. power_A ==0)then
double precision :: F_integral
overlap_x = P_center * F_integral(0,p) * factor
dx = 0.d0
return
endif
double precision :: pouet_timy
pouet_timy = dsqrt(lower_exp_val/p)

2
src/NEEDED_MODULES
View File

@ -1 +1 @@
AOs BiInts Bitmask Dets Electrons Ezfio_files Full_CI Generators_full Hartree_Fock MOGuess MonoInts MOs Nuclei Output Selectors_full Utils Molden FCIdump Generators_CAS CAS_SD_selected DDCI_selected
AOs BiInts Bitmask Dets Electrons Ezfio_files Full_CI Generators_full Hartree_Fock MOGuess MonoInts MOs Nuclei Output Selectors_full Utils Molden FCIdump Generators_CAS CAS_SD_selected

16
src/Nuclei/nuclei.irp.f
View File

@ -168,6 +168,22 @@ END_PROVIDER
END_PROVIDER
BEGIN_PROVIDER [ double precision, positive_charge_barycentre,(3)]
implicit none
BEGIN_DOC
! Centroid of the positive charges
END_DOC
integer :: l
positive_charge_barycentre(1) = 0.d0
positive_charge_barycentre(2) = 0.d0
positive_charge_barycentre(3) = 0.d0
do l = 1, nucl_num
positive_charge_barycentre(1) += nucl_charge(l) * nucl_coord(l,1)
positive_charge_barycentre(2) += nucl_charge(l) * nucl_coord(l,2)
positive_charge_barycentre(3) += nucl_charge(l) * nucl_coord(l,3)
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, nuclear_repulsion ]
implicit none
BEGIN_DOC

2
src/Perturbation/NEEDED_MODULES
View File

@ -1,2 +1,2 @@
AOs BiInts Bitmask Dets Electrons Ezfio_files Hartree_Fock MOGuess MonoInts MOs Nuclei Output Utils
AOs BiInts Bitmask Dets Electrons Ezfio_files Hartree_Fock MOGuess MonoInts MOs Nuclei Output Properties Utils

72
src/Perturbation/delta_rho_perturbation.irp.f Normal file
View File

@ -0,0 +1,72 @@
subroutine pt2_delta_rho_one_point(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st)
use bitmasks
implicit none
integer, intent(in) :: Nint,ndet,n_st
integer(bit_kind), intent(in) :: det_pert(Nint,2)
double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag(N_st)
double precision :: i_O1_psi_array(N_st)
double precision :: i_H_psi_array(N_st)
BEGIN_DOC
! compute the perturbatibe contribution to the Integrated Spin density at z = z_one point of one determinant
!
! for the various n_st states, at various level of theory.
!
! c_pert(i) = <psi(i)|H|det_pert>/(<psi(i)|H|psi(i)> - <det_pert|H|det_pert>)
!
! e_2_pert(i) = c_pert(i) * <det_pert|O|psi(i)>
!
! H_pert_diag(i) = c_pert(i)^2 * <det_pert|O|det_pert>
!
! To get the contribution of the first order :
!
! <O_1> = sum(over i) e_2_pert(i)
!
! To get the contribution of the diagonal elements of the second order :
!
! [ <O_0> + <O_1> + sum(over i) H_pert_diag(i) ] / [1. + sum(over i) c_pert(i) **2]
!
END_DOC
integer :: i,j
double precision :: diag_H_mat_elem,diag_o1_mat_elem_alpha_beta
integer :: exc(0:2,2,2)
integer :: degree
double precision :: phase,delta_e,h,oii,diag_o1_mat_elem
integer :: h1,h2,p1,p2,s1,s2
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
! call get_excitation_degree(HF_bitmask,det_pert,degree,N_int)
! if(degree.gt.degree_max_generators+1)then
! H_pert_diag = 0.d0
! e_2_pert = 0.d0
! c_pert = 0.d0
! return
! endif
call i_O1_psi_alpha_beta(mo_integrated_delta_rho_one_point,det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_O1_psi_array)
call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
h = diag_H_mat_elem(det_pert,Nint)
oii = diag_O1_mat_elem_alpha_beta(mo_integrated_delta_rho_one_point,det_pert,N_int)
do i =1,N_st
if(CI_electronic_energy(i)>h.and.CI_electronic_energy(i).ne.0.d0)then
c_pert(i) = -1.d0
e_2_pert(i) = selection_criterion*selection_criterion_factor*2.d0
else if (dabs(CI_electronic_energy(i) - h) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - h)
e_2_pert(i) = c_pert(i) * (i_O1_psi_array(i)+i_O1_psi_array(i) ) + c_pert(i) * c_pert(i) * oii
H_pert_diag(i) = c_pert(i) * (i_O1_psi_array(i)+i_O1_psi_array(i) )
else
c_pert(i) = -1.d0
e_2_pert(i) = -dabs(i_H_psi_array(i))
H_pert_diag(i) = c_pert(i) * i_O1_psi_array(i)
endif
enddo
end

69
src/Perturbation/dipole_moment.irp.f Normal file
View File

@ -0,0 +1,69 @@
subroutine pt2_dipole_moment_z(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st)
use bitmasks
implicit none
integer, intent(in) :: Nint,ndet,n_st
integer(bit_kind), intent(in) :: det_pert(Nint,2)
double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag(N_st)
double precision :: i_O1_psi_array(N_st)
double precision :: i_H_psi_array(N_st)
BEGIN_DOC
! compute the perturbatibe contribution to the dipole moment of one determinant
!
! for the various n_st states, at various level of theory.
!
! c_pert(i) = <psi(i)|H|det_pert>/(<psi(i)|H|psi(i)> - <det_pert|H|det_pert>)
!
! e_2_pert(i) = c_pert(i) * <det_pert|Z|psi(i)>
!
! H_pert_diag(i) = c_pert(i)^2 * <det_pert|Z|det_pert>
!
! To get the contribution of the first order :
!
! <Z_1> = sum(over i) e_2_pert(i)
!
! To get the contribution of the diagonal elements of the second order :
!
! [ <Z_0> + <Z_1> + sum(over i) H_pert_diag(i) ] / [1. + sum(over i) c_pert(i) **2]
!
END_DOC
integer :: i,j
double precision :: diag_H_mat_elem
integer :: exc(0:2,2,2)
integer :: degree
double precision :: phase,delta_e,h,oii,diag_o1_mat_elem
integer :: h1,h2,p1,p2,s1,s2
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
! call get_excitation_degree(HF_bitmask,det_pert,degree,N_int)
! if(degree.gt.degree_max_generators+1)then
! H_pert_diag = 0.d0
! e_2_pert = 0.d0
! c_pert = 0.d0
! return
! endif
call i_O1_psi(mo_dipole_z,det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_O1_psi_array)
call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
h = diag_H_mat_elem(det_pert,Nint)
oii = diag_O1_mat_elem(mo_dipole_z,det_pert,N_int)
do i =1,N_st
if(CI_electronic_energy(i)>h.and.CI_electronic_energy(i).ne.0.d0)then
c_pert(i) = -1.d0
e_2_pert(i) = selection_criterion*selection_criterion_factor*2.d0
else if (dabs(CI_electronic_energy(i) - h) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - h)
e_2_pert(i) = c_pert(i) * (i_O1_psi_array(i)+i_O1_psi_array(i) )
H_pert_diag(i) = e_2_pert(i) + c_pert(i) * c_pert(i) * oii
else
c_pert(i) = -1.d0
e_2_pert(i) = -dabs(i_H_psi_array(i))
H_pert_diag(i) = c_pert(i) * i_O1_psi_array(i)
endif
enddo
end

4
src/Perturbation/exc_max.irp.f Normal file
View File

@ -0,0 +1,4 @@
BEGIN_PROVIDER [integer, max_exc_pert]
implicit none
max_exc_pert = 0
END_PROVIDER

49
src/Perturbation/pert_single.irp.f Normal file
View File

@ -0,0 +1,49 @@
subroutine pt2_h_core(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st)
use bitmasks
implicit none
integer, intent(in) :: Nint,ndet,N_st
integer(bit_kind), intent(in) :: det_pert(Nint,2)
double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st)
double precision :: i_H_psi_array(N_st)
BEGIN_DOC
! compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
!
! for the various N_st states.
!
! c_pert(i) = <psi(i)|H|det_pert>/( E(i) - <det_pert|H|det_pert> )
!
! e_2_pert(i) = <psi(i)|H|det_pert>^2/( E(i) - <det_pert|H|det_pert> )
!
END_DOC
integer :: i,j
double precision :: diag_H_mat_elem, h
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
integer :: exc(0:2,2,2)
integer :: degree
double precision :: phase
call get_excitation(ref_bitmask,det_pert,exc,degree,phase,N_int)
h = diag_H_mat_elem(det_pert,N_int)
print*,'delta E = ',h-ref_bitmask_energy
if(h<ref_bitmask_energy)then
c_pert = 0.d0
e_2_pert = 0.d0
H_pert_diag = 0.d0
return
endif
if(degree>1)then
c_pert = 0.d0
e_2_pert = 0.d0
H_pert_diag = 0.d0
return
endif
integer :: h1,p1,h2,p2,s1,s2
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
c_pert = phase * mo_mono_elec_integral(h1,p1)
e_2_pert = -dabs(mo_mono_elec_integral(h1,p1)+1.d0)
end

52
src/Perturbation/perturbation_template.f
View File

@ -35,6 +35,58 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c
cycle
endif
integer :: degree
call get_excitation_degree(HF_bitmask,buffer(1,1,i),degree,N_int)
call pt2_$PERT(buffer(1,1,i), &
c_pert,e_2_pert,H_pert_diag,Nint,N_det_selectors,n_st)
do k = 1,N_st
e_2_pert_buffer(k,i) = e_2_pert(k)
coef_pert_buffer(k,i) = c_pert(k)
sum_norm_pert(k) += c_pert(k) * c_pert(k)
sum_e_2_pert(k) += e_2_pert(k)
sum_H_pert_diag(k) += H_pert_diag(k)
enddo
enddo
end
subroutine perturb_buffer_by_mono_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint)
implicit none
BEGIN_DOC
! Applly pertubration ``$PERT`` to the buffer of determinants generated in the H_apply
! routine.
END_DOC
integer, intent(in) :: Nint, N_st, buffer_size, i_generator
integer(bit_kind), intent(in) :: buffer(Nint,2,buffer_size)
double precision, intent(inout) :: sum_norm_pert(N_st),sum_e_2_pert(N_st)
double precision, intent(inout) :: coef_pert_buffer(N_st,buffer_size),e_2_pert_buffer(N_st,buffer_size),sum_H_pert_diag(N_st)
double precision :: c_pert(N_st), e_2_pert(N_st), H_pert_diag(N_st)
integer :: i,k, c_ref
integer, external :: connected_to_ref_by_mono
logical, external :: is_in_wavefunction
ASSERT (Nint > 0)
ASSERT (Nint == N_int)
ASSERT (buffer_size >= 0)
ASSERT (minval(sum_norm_pert) >= 0.d0)
ASSERT (N_st > 0)
do i = 1,buffer_size
c_ref = connected_to_ref_by_mono(buffer(1,1,i),psi_generators,Nint,i_generator,N_det)
if (c_ref /= 0) then
cycle
endif
if (is_in_wavefunction(buffer(1,1,i),Nint,N_det)) then
cycle
endif
integer :: degree
call pt2_$PERT(buffer(1,1,i), &
c_pert,e_2_pert,H_pert_diag,Nint,N_det_selectors,n_st)

4
src/Perturbation/selection.irp.f
View File

@ -36,7 +36,7 @@ subroutine fill_H_apply_buffer_selection(n_selected,det_buffer,e_2_pert_buffer,c
is_selected = .False.
do j=1,N_st
s = -e_2_pert_buffer(j,i)
s = dabs(e_2_pert_buffer(j,i))
is_selected = s > selection_criterion*selection_criterion_factor .or. is_selected
select_max_out = max(select_max_out,s)
enddo
@ -78,7 +78,7 @@ end
BEGIN_DOC
! Threshold to select determinants. Set by selection routines.
END_DOC
selection_criterion = .1d0
selection_criterion = 0.1d0
selection_criterion_factor = 0.01d0
selection_criterion_min = selection_criterion

1
src/Primitive_basis/NEEDED_MODULES
View File

@ -1 +0,0 @@
AOs Electrons Ezfio_files MOs Nuclei Output Utils

18
src/Primitive_basis/functions.irp.f
View File

@ -1,18 +0,0 @@
double precision function mo_r(i,r)
implicit none
integer :: i
double precision :: r(3)
integer :: j
double precision :: x_center,y_center,z_center,r_center
mo_r = 0.d0
do j = 1, prim_num
if(dabs(prim_mo_coef(j,i)).lt.1.d-10)cycle
x_center = r(1) - nucl_coord(prim_nucl(j),1)
y_center = r(2) - nucl_coord(prim_nucl(j),2)
z_center = r(3) - nucl_coord(prim_nucl(j),3)
r_center = x_center*x_center + y_center*y_center + z_center*z_center
mo_r = mo_r + prim_mo_coef(j,i) * (x_center ** (prim_power(j,1)) * y_center ** (prim_power(j,2)) * z_center ** (prim_power(j,3)) &
* dexp(-prim_expo(j) * r_center))
enddo
end

175
src/Primitive_basis/prim.irp.f
View File

@ -1,175 +0,0 @@
BEGIN_PROVIDER [integer, prim_num]
implicit none
BEGIN_DOC
! Number of uniq primitive basis function
END_DOC
PROVIDE ezfio_filename
call ezfio_get_primitive_basis_prim_num(prim_num)
END_PROVIDER
BEGIN_PROVIDER [integer, prim_nucl, (prim_num)]
implicit none
BEGIN_DOC
! Array of the nuclei on which are attached the primitives
END_DOC
PROVIDE ezfio_filename
integer, allocatable :: tmp(:)
allocate(tmp(prim_num))
call ezfio_get_primitive_basis_prim_nucl(tmp)
prim_nucl = tmp
END_PROVIDER
BEGIN_PROVIDER [integer, prim_power, (prim_num,3)]
&BEGIN_PROVIDER [integer, prim_l, (prim_num)]
implicit none
BEGIN_DOC
! Array of the power of the primitives
! Array of the L of the primitive
END_DOC
PROVIDE ezfio_filename
integer :: i,j
integer, allocatable :: tmp(:,:)
allocate(tmp(3,prim_num))
call ezfio_get_primitive_basis_prim_power(tmp)
do i = 1, prim_num
do j = 1,3
prim_power(i,j) = tmp(j,i)
enddo
enddo
do i = 1, prim_num
prim_l(i) = prim_power(i,1) + prim_power(i,2) + prim_power(i,3)
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, prim_expo, (prim_num)]
implicit none
BEGIN_DOC
! Array of the exponents of the primitives
END_DOC
PROVIDE ezfio_filename
double precision, allocatable :: tmp(:)
allocate(tmp(prim_num))
call ezfio_get_primitive_basis_prim_expo(tmp)
prim_expo = tmp
END_PROVIDER
BEGIN_PROVIDER [double precision, prim_overlap, (prim_num,prim_num)]
implicit none
BEGIN_DOC
! Array of the overlap between the primitives
END_DOC
integer :: i,j
prim_overlap = 0.d0
double precision :: overlap, overlap_x, overlap_y, overlap_z
double precision :: alpha, beta, c
double precision :: A_center(3), B_center(3)
integer :: power_A(3), power_B(3),dim1
dim1=100
do j = 1, prim_num
A_center(1) = nucl_coord( prim_nucl(j), 1 )
A_center(2) = nucl_coord( prim_nucl(j), 2 )
A_center(3) = nucl_coord( prim_nucl(j), 3 )
power_A(1) = prim_power( j, 1 )
power_A(2) = prim_power( j, 2 )
power_A(3) = prim_power( j, 3 )
alpha = prim_expo(j)
do i = 1, prim_num
B_center(1) = nucl_coord( prim_nucl(i), 1 )
B_center(2) = nucl_coord( prim_nucl(i), 2 )
B_center(3) = nucl_coord( prim_nucl(i), 3 )
power_B(1) = prim_power( i, 1 )
power_B(2) = prim_power( i, 2 )
power_B(3) = prim_power( i, 3 )
beta = prim_expo(i)
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1)
prim_overlap(i,j) = overlap
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, prim_ao_coef, (prim_num,ao_num)]
BEGIN_DOC
! Developement of the primitive on the AO basis
! prim_ao_coef(i,j) = <prim_i|AO_j>
END_DOC
implicit none
integer :: i,j,k,l
prim_ao_coef = 0.d0
do i = 1, prim_num
do j = 1, ao_num
if(ao_nucl(j).ne.prim_nucl(i))cycle
if(ao_l(j) == prim_l(i))then
if(ao_power(j,1) == prim_power(i,1) .and. &
ao_power(j,2) == prim_power(i,2) .and. &
ao_power(j,3) == prim_power(i,3) )then
! primitive of the same L than the AO
do k = 1, ao_prim_num(j)
if(ao_expo(j,k) == prim_expo(i))then
prim_ao_coef(i,j) += ao_coef(j,k)
endif
enddo
endif
endif
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, prim_ao_overlap, (ao_num,ao_num)]
implicit none
BEGIN_dOC
! Array of the overlap of the AO basis, explicitely calculated on the primitive basis
END_DOC
integer :: i,j,k,l
do i= 1, ao_num
do j = 1,ao_num
prim_ao_overlap(i,j) = 0.d0
do k = 1, prim_num
do l = 1, prim_num
prim_ao_overlap(i,j) += prim_ao_coef(k,i) * prim_overlap(k,l) * prim_ao_coef(l,j)
enddo
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, prim_mo_coef, (prim_num,mo_tot_num)]
BEGIN_DOC
! Developement of the primitive on the MO basis
! prim_ao_coef(i,j) = <prim(i)|ao(j)>
END_DOC
integer :: i,j,k,l
double precision :: coef
prim_mo_coef = 0.d0
do i = 1, mo_tot_num
do k = 1, ao_num
coef = mo_coef(k,i)
do l = 1, prim_num
prim_mo_coef(l,i) += coef * prim_ao_coef(l,k)
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, prim_mo_overlap, (mo_tot_num,mo_tot_num)]
implicit none
BEGIN_dOC
! Array of the overlap of the AO basis, explicitely calculated on the primitive basis
END_DOC
integer :: i,j,k,l
do i= 1, mo_tot_num
do j = 1,mo_tot_num
prim_mo_overlap(i,j) = 0.d0
do k = 1, prim_num
do l = 1, prim_num
prim_mo_overlap(i,j) += prim_mo_coef(k,i) * prim_overlap(k,l) * prim_mo_coef(l,j)
enddo
enddo
enddo
enddo
END_PROVIDER

5
src/Primitive_basis/primitive_basis.ezfio_config
View File

@ -1,5 +0,0 @@
primitive_basis
prim_num integer
prim_power integer (3,primitive_basis_prim_num)
prim_expo double precision (primitive_basis_prim_num)
prim_nucl integer (primitive_basis_prim_num)

8
src/Properties/Makefile Normal file
View File

@ -0,0 +1,8 @@
default: all
# Define here all new external source files and objects.Don't forget to prefix the
# object files with IRPF90_temp/
SRC=
OBJ=
include $(QPACKAGE_ROOT)/src/Makefile.common

1
src/Properties/NEEDED_MODULES Normal file
View File

@ -0,0 +1 @@
AOs BiInts Bitmask Dets Electrons Ezfio_files MonoInts MOs Nuclei Output Utils

4
src/Properties/README.rst Normal file
View File

@ -0,0 +1,4 @@
=================
Properties Module
=================

25
src/Properties/average.irp.f Normal file
View File

@ -0,0 +1,25 @@
subroutine get_average(array,density,average)
implicit none
double precision, intent(in) :: array(mo_tot_num_align,mo_tot_num)
double precision, intent(in) :: density(mo_tot_num_align,mo_tot_num)
double precision, intent(out):: average
integer :: i,j
BEGIN_DOC
! computes the average value of a pure MONO ELECTRONIC OPERATOR
! whom integrals on the MO basis are stored in "array"
! and with the density is stored in "density"
END_DOC
average = 0.d0
!$OMP PARALLEL DO DEFAULT(none) &
!$OMP PRIVATE(i,j) &
!$OMP SHARED(mo_tot_num,array,density) &
!$OMP REDUCTION(+:average)
do i = 1, mo_tot_num
do j = 1, mo_tot_num
average += density(j,i) * array(j,i)
enddo
enddo
!$OMP END PARALLEL DO
end

224
src/Properties/delta_rho.irp.f Normal file
View File

@ -0,0 +1,224 @@
BEGIN_PROVIDER [integer, N_z_pts]
&BEGIN_PROVIDER [double precision, z_min]
&BEGIN_PROVIDER [double precision, z_max]
&BEGIN_PROVIDER [double precision, delta_z]
implicit none
z_min = -20.d0
z_max = 20.d0
delta_z = 0.1d0
N_z_pts = (z_max - z_min)/delta_z
print*,'N_z_pts = ',N_z_pts
END_PROVIDER
BEGIN_PROVIDER [double precision, integrated_delta_rho_all_points, (N_z_pts)]
BEGIN_DOC
!
! integrated_rho(alpha,z) - integrated_rho(beta,z) for all the z points
! chosen
!
END_DOC
implicit none
integer :: i,j,k,l,i_z,h
double precision :: z,function_integrated_delta_rho,c_k,c_j,n_i_h,accu
integrated_delta_rho_all_points = 0.d0
!$OMP PARALLEL DO DEFAULT(none) &
!$OMP PRIVATE(i,h,j,k,c_j,c_k,n_i_h,i_z) &
!$OMP SHARED(mo_tot_num,ao_num,mo_coef, &
!$OMP ao_integrated_delta_rho_all_points,one_body_spin_density_mo,integrated_delta_rho_all_points,N_z_pts)
do i_z = 1, N_z_pts
do i = 1, mo_tot_num
do h = 1, mo_tot_num
n_i_h = one_body_spin_density_mo(i,h)
if(dabs(n_i_h).lt.1.d-10)cycle
do j = 1, ao_num
c_j = mo_coef(j,i) ! coefficient of the ith MO on the jth AO
do k = 1, ao_num
c_k = mo_coef(k,h) ! coefficient of the hth MO on the kth AO
integrated_delta_rho_all_points(i_z) += c_k * c_j * n_i_h * ao_integrated_delta_rho_all_points(j,k,i_z)
enddo
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
z = z_min
accu = 0.d0
do i = 1, N_z_pts
accu += integrated_delta_rho_all_points(i)
write(i_unit_integrated_delta_rho,*)z,integrated_delta_rho_all_points(i),accu
z += delta_z
enddo
print*,'sum of integrated_delta_rho = ',accu
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_integrated_delta_rho_all_points, (ao_num_align, ao_num, N_z_pts)]
BEGIN_DOC
! array of the overlap in x,y between the AO function and integrated between [z,z+dz] in the z axis
! for all the z points that are given (N_z_pts)
END_DOC
implicit none
integer :: i,j,n,l
double precision :: f,accu
integer :: dim1
double precision :: overlap, overlap_x, overlap_y, overlap_z
double precision :: alpha, beta, c
double precision :: A_center(3), B_center(3)
integer :: power_A(3), power_B(3)
integer :: i_z
double precision :: z,SABpartial,accu_x,accu_y,accu_z
dim1=100
z = z_min
do i_z = 1, N_z_pts
!$OMP PARALLEL DO DEFAULT(none) &
!$OMP PRIVATE(i,j,n,l,A_center,power_A,B_center,power_B,accu_z, &
!$OMP overlap_x,overlap_y,overlap_z,overlap,c,alpha,beta) &
!$OMP SHARED(ao_num,nucl_coord,ao_nucl,ao_power,ao_prim_num,ao_expo_transp,ao_coef_transp, &
!$OMP ao_integrated_delta_rho_all_points,N_z_pts,dim1,i_z,z,delta_z)
do j=1,ao_num
A_center(1) = nucl_coord( ao_nucl(j), 1 )
A_center(2) = nucl_coord( ao_nucl(j), 2 )
A_center(3) = nucl_coord( ao_nucl(j), 3 )
power_A(1) = ao_power( j, 1 )
power_A(2) = ao_power( j, 2 )
power_A(3) = ao_power( j, 3 )
do i= 1,ao_num
B_center(1) = nucl_coord( ao_nucl(i), 1 )
B_center(2) = nucl_coord( ao_nucl(i), 2 )
B_center(3) = nucl_coord( ao_nucl(i), 3 )
power_B(1) = ao_power( i, 1 )
power_B(2) = ao_power( i, 2 )
power_B(3) = ao_power( i, 3 )
accu_z = 0.d0
do n = 1,ao_prim_num(j)
alpha = ao_expo_transp(n,j)
do l = 1, ao_prim_num(i)
beta = ao_expo_transp(l,i)
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1)
c = ao_coef_transp(n,j) * ao_coef_transp(l,i)
accu_z += c* overlap_x * overlap_y * SABpartial(z,z+delta_z,A_center,B_center,power_A,power_B,alpha,beta)
enddo
enddo
ao_integrated_delta_rho_all_points(i,j,i_z) = accu_z
enddo
enddo
!$OMP END PARALLEL DO
z += delta_z
enddo
END_PROVIDER
BEGIN_PROVIDER [integer, i_unit_integrated_delta_rho]
implicit none
BEGIN_DOC
! fortran unit for the writing of the integrated delta_rho
END_DOC
integer :: getUnitAndOpen
character*(128) :: output_i_unit_integrated_delta_rho
output_i_unit_integrated_delta_rho=trim(ezfio_filename)//'/properties/delta_rho'
i_unit_integrated_delta_rho= getUnitAndOpen(output_i_unit_integrated_delta_rho,'w')
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_integrated_delta_rho_one_point, (ao_num_align, ao_num )]
BEGIN_DOC
! array of the overlap in x,y between the AO function and integrated between [z,z+dz] in the z axis
! for one specific z point
END_DOC
implicit none
integer :: i,j,n,l
double precision :: f
integer :: dim1
double precision :: overlap, overlap_x, overlap_y, overlap_z
double precision :: alpha, beta, c
double precision :: A_center(3), B_center(3)
integer :: power_A(3), power_B(3)
integer :: i_z
double precision :: z,SABpartial,accu_z
dim1=100
z = z_one_point
!$OMP PARALLEL DO DEFAULT(none) &
!$OMP PRIVATE(i,j,n,l,A_center,power_A,B_center,power_B,accu_z, &
!$OMP overlap_x,overlap_y,overlap_z,overlap,c,alpha,beta) &
!$OMP SHARED(ao_num,nucl_coord,ao_nucl,ao_power,ao_prim_num,ao_expo_transp,ao_coef_transp, &
!$OMP ao_integrated_delta_rho_one_point,dim1,z,delta_z)
do j=1,ao_num
A_center(1) = nucl_coord( ao_nucl(j), 1 )
A_center(2) = nucl_coord( ao_nucl(j), 2 )
A_center(3) = nucl_coord( ao_nucl(j), 3 )
power_A(1) = ao_power( j, 1 )
power_A(2) = ao_power( j, 2 )
power_A(3) = ao_power( j, 3 )
do i= 1,ao_num
B_center(1) = nucl_coord( ao_nucl(i), 1 )
B_center(2) = nucl_coord( ao_nucl(i), 2 )
B_center(3) = nucl_coord( ao_nucl(i), 3 )
power_B(1) = ao_power( i, 1 )
power_B(2) = ao_power( i, 2 )
power_B(3) = ao_power( i, 3 )
accu_z = 0.d0
do n = 1,ao_prim_num(j)
alpha = ao_expo_transp(n,j)
do l = 1, ao_prim_num(i)
beta = ao_expo_transp(l,i)
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1)
c = ao_coef_transp(n,j) * ao_coef_transp(l,i)
accu_z += c* overlap_x * overlap_y * SABpartial(z,z+delta_z,A_center,B_center,power_A,power_B,alpha,beta)
enddo
enddo
ao_integrated_delta_rho_one_point(i,j) = accu_z
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER
BEGIN_PROVIDER [double precision, mo_integrated_delta_rho_one_point, (mo_tot_num_align,mo_tot_num)]
BEGIN_DOC
!
! array of the integrals needed of integrated_rho(alpha,z) - integrated_rho(beta,z) for z = z_one_point
! on the MO basis
!
END_DOC
implicit none
integer :: i,j,k,l,i_z,h
double precision :: z,function_integrated_delta_rho,c_k,c_j
mo_integrated_delta_rho_one_point = 0.d0
!$OMP PARALLEL DO DEFAULT(none) &
!$OMP PRIVATE(i,j,h,k,c_j,c_k) &
!$OMP SHARED(mo_tot_num,ao_num,mo_coef, &
!$OMP mo_integrated_delta_rho_one_point, ao_integrated_delta_rho_one_point)
do i = 1, mo_tot_num
do h = 1, mo_tot_num
do j = 1, ao_num
c_j = mo_coef(j,i) ! coefficient of the jth AO on the ith MO
do k = 1, ao_num
c_k = mo_coef(k,h) ! coefficient of the kth AO on the hth MO
mo_integrated_delta_rho_one_point(i,h) += c_k * c_j * ao_integrated_delta_rho_one_point(j,k)
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER
BEGIN_PROVIDER [ double precision, integrated_delta_rho_one_point]
implicit none
BEGIN_DOC
!
! integral (x,y) and (z,z+delta_z) of rho(alpha) - rho(beta)
! on the MO basis
!
END_DOC
double precision :: average
call get_average(mo_integrated_delta_rho_one_point,one_body_spin_density_mo,average)
integrated_delta_rho_one_point = average
END_PROVIDER

289
src/Properties/need.irp.f Normal file
View File

@ -0,0 +1,289 @@
double precision function SABpartial(zA,zB,A,B,nA,nB,gamA,gamB)
implicit double precision(a-h,o-z)
dimension nA(3),nB(3)
dimension A(3),B(3)
gamtot=gamA+gamB
SABpartial=1.d0
l=3
u=gamA/gamtot*A(l)+gamB/gamtot*B(l)
arg=gamtot*u**2-gamA*A(l)**2-gamB*B(l)**2
alpha=dexp(arg)
&/gamtot**((1.d0+dfloat(nA(l))+dfloat(nB(l)))/2.d0)
wA=dsqrt(gamtot)*(u-A(l))
wB=dsqrt(gamtot)*(u-B(l))
boundA=dsqrt(gamtot)*(zA-u)
boundB=dsqrt(gamtot)*(zB-u)
accu=0.d0
do n=0,nA(l)
do m=0,nB(l)
integ=nA(l)+nB(l)-n-m
accu=accu
& +wA**n*wB**m*binom(nA(l),n)*binom(nB(l),m)
& *(rinteg(integ,boundB)-rinteg(integ,boundA))
enddo
enddo
SABpartial=SABpartial*accu*alpha
end
double precision function rintgauss(n)
implicit double precision(a-h,o-z)
rintgauss=dsqrt(dacos(-1.d0))
if(n.eq.0)return
if(n.eq.1)then
rintgauss=0.d0
return
endif
if(iand(n,1).eq.1)then
rintgauss=0.d0
return
endif
rintgauss=rintgauss/2.d0**(n/2)
rintgauss=rintgauss*ddfact2(n-1)
end
double precision function rinteg(n,u)
implicit double precision(a-h,o-z)
include 'constants.F'
! pi=dacos(-1.d0)
ichange=1
factor=1.d0
if(u.lt.0.d0)then
u=-u
factor=(-1.d0)**(n+1)
ichange=-1
endif
if(iand(n,1).eq.0)then
rinteg=0.d0
do l=0,n-2,2
prod=b_coef(l,u)
do k=l+2,n-2,2
prod=prod*a_coef(k)
enddo
rinteg=rinteg+prod
enddo
prod=dsqrt(pi)/2.d0*erf0(u)
do k=0,n-2,2
prod=prod*a_coef(k)
enddo
rinteg=rinteg+prod
endif
if(iand(n,1).eq.1)then
rinteg=0.d0
do l=1,n-2,2
prod=b_coef(l,u)
do k=l+2,n-2,2
prod=prod*a_coef(k)
enddo
rinteg=rinteg+prod
enddo
prod=0.5d0*(1.d0-dexp(-u**2))
do k=1,n-2,2
prod=prod*a_coef(k)
enddo
rinteg=rinteg+prod
endif
rinteg=rinteg*factor
if(ichange.eq.-1)u=-u
end
!<function type="double precision function" name="erf0">
! <arg name="x"
! doc ="" />
!
! <doc>
!
! </doc>
!
! <fortran>
double precision function erf0(x)
implicit double precision (a-h,o-z)
if(x.lt.0.d0)then
erf0=-gammp(0.5d0,x**2)
else
erf0=gammp(0.5d0,x**2)
endif
end
! </fortran>
!</function>
!<function type="double precision function" name="gammp">
! <arg name="a"
! doc ="" />
! <arg name="x"
! doc ="" />
!
! <doc>
!
! </doc>
!
! <calls>
! gcf
! gser
! </calls>
!
! <fortran>
double precision function gammp(a,x)
implicit double precision (a-h,o-z)
if(x.lt.0..or.a.le.0.)pause
if(x.lt.a+1.)then
call gser(gammp,a,x,gln)
else
call gcf(gammcf,a,x,gln)
gammp=1.-gammcf
endif
return
end
! </fortran>
!</function>
!<function type="subroutine" name="gser">
! <arg name="gamser"
! doc ="" />
! <arg name="a"
! doc ="" />
! <arg name="x"
! doc ="" />
! <arg name="gln"
! doc ="" />
!
! <doc>
!
! </doc>
!
! <calledBy>
! gammp
! </calledBy>
!
! <fortran>
subroutine gser(gamser,a,x,gln)
implicit double precision (a-h,o-z)
parameter (itmax=100,eps=3.e-7)
gln=gammln(a)
if(x.le.0.)then
if(x.lt.0.)pause
gamser=0.
return
endif
ap=a
sum=1./a
del=sum
do 11 n=1,itmax
ap=ap+1.
del=del*x/ap
sum=sum+del
if(abs(del).lt.abs(sum)*eps)go to 1
11 continue
pause 'a too large, itmax too small'
1 gamser=sum*exp(-x+a*log(x)-gln)
return
end
! </fortran>
!</function>
!<function type="subroutine" name="gcf">
! <arg name="gammcf"
! doc ="" />
! <arg name="a"
! doc ="" />
! <arg name="x"
! doc ="" />
! <arg name="gln"
! doc ="" />
!
! <doc>
!
! </doc>
!
! <calledBy>
! gammp
! </calledBy>
!
! <fortran>
subroutine gcf(gammcf,a,x,gln)
implicit double precision (a-h,o-z)
parameter (itmax=100,eps=3.e-7)
gln=gammln(a)
gold=0.
a0=1.
a1=x
b0=0.
b1=1.
fac=1.
do 11 n=1,itmax
an=float(n)
ana=an-a
a0=(a1+a0*ana)*fac
b0=(b1+b0*ana)*fac
anf=an*fac
a1=x*a0+anf*a1
b1=x*b0+anf*b1
if(a1.ne.0.)then
fac=1./a1
g=b1*fac
if(abs((g-gold)/g).lt.eps)go to 1
gold=g
endif
11 continue
pause 'a too large, itmax too small'
1 gammcf=exp(-x+a*log(x)-gln)*g
return
end
! </fortran>
!</function>
double precision function ddfact2(n)
implicit double precision(a-h,o-z)
if(iand(n,1).eq.0)stop 'error in ddfact2'
ddfact2=1.d0
do i=1,n,2
ddfact2=ddfact2*dfloat(i)
enddo
end
double precision function a_coef(n)
implicit double precision(a-h,o-z)
a_coef=dfloat(n+1)/2.d0
end
double precision function b_coef(n,u)
implicit double precision(a-h,o-z)
b_coef=-0.5d0*u**(n+1)*dexp(-u**2)
end
!<function type="double precision function" name="gammln">
! <arg name="xx"
! doc ="" />
!
! <doc>
!
! </doc>
!
! <fortran>
double precision function gammln(xx)
implicit double precision (a-h,o-z)
real*8 cof(6),stp,half,one,fpf,x,tmp,ser
data cof,stp/76.18009173d0,-86.50532033d0,24.01409822d0,
* -1.231739516d0,.120858003d-2,-.536382d-5,2.50662827465d0/
data half,one,fpf/0.5d0,1.0d0,5.5d0/
x=xx-one
tmp=x+fpf
tmp=(x+half)*log(tmp)-tmp
ser=one
do 11 j=1,6
x=x+one
ser=ser+cof(j)/x
11 continue
gammln=tmp+log(stp*ser)
return
end
! </fortran>
!</function>

13
src/Properties/options.irp.f Normal file
View File

@ -0,0 +1,13 @@
BEGIN_SHELL [ /usr/bin/python ]
from ezfio_with_default import EZFIO_Provider
T = EZFIO_Provider()
T.set_type ( "double precision" )
T.set_name ( "z_one_point" )
T.set_doc ( "z point on which the integrated delta rho is calculated" )
T.set_ezfio_dir ( "properties" )
T.set_ezfio_name( "z_one_point" )
T.set_output ( "output_full_ci" )
print T
END_SHELL

2
src/Properties/properties.ezfio_config Normal file
View File

@ -0,0 +1,2 @@
properties
z_one_point double precision

51
src/Properties/properties.irp.f Normal file
View File

@ -0,0 +1,51 @@
BEGIN_PROVIDER [double precision, average_position,(3)]
implicit none
BEGIN_DOC
! average_position(1) = <psi_det|X|psi_det>
! average_position(2) = <psi_det|Y|psi_det>
! average_position(3) = <psi_det|Z|psi_det>
END_DOC
integer :: i,j
average_position = 0.d0
!$OMP PARALLEL DO DEFAULT(none) &
!$OMP PRIVATE(i,j) &
!$OMP SHARED(mo_tot_num,mo_dipole_x,mo_dipole_y,mo_dipole_z,one_body_dm_mo) &
!$OMP REDUCTION(+:average_position)
do i = 1, mo_tot_num
do j = 1, mo_tot_num
average_position(1) += one_body_dm_mo(j,i) * mo_dipole_x(j,i)
average_position(2) += one_body_dm_mo(j,i) * mo_dipole_y(j,i)
average_position(3) += one_body_dm_mo(j,i) * mo_dipole_z(j,i)
enddo
enddo
!$OMP END PARALLEL DO
!call test_average_value(mo_dipole_z,average_position(3))
END_PROVIDER
BEGIN_PROVIDER [double precision, average_spread,(3)]
implicit none
BEGIN_DOC
! average_spread(1) = <psi_det|X^2|psi_det>
! average_spread(2) = <psi_det|Y^2|psi_det>
! average_spread(3) = <psi_det|Z^2|psi_det>
END_DOC
integer :: i,j
average_spread = 0.d0
!$OMP PARALLEL DO DEFAULT(none) &
!$OMP PRIVATE(i,j) &
!$OMP SHARED(mo_tot_num,mo_spread_x,mo_spread_y,mo_spread_z,one_body_dm_mo) &
!$OMP REDUCTION(+:average_spread)
do i = 1, mo_tot_num
do j = 1, mo_tot_num
average_spread(1) += one_body_dm_mo(j,i) * mo_spread_x(j,i)
average_spread(2) += one_body_dm_mo(j,i) * mo_spread_y(j,i)
average_spread(3) += one_body_dm_mo(j,i) * mo_spread_z(j,i)
enddo
enddo
!$OMP END PARALLEL DO
!call test_average_value(mo_spread_z,average_spread(3))
END_PROVIDER

89
src/Properties/routines_test.irp.f Normal file
View File

@ -0,0 +1,89 @@
subroutine test_average_value(array,value)
implicit none
double precision, intent(in) :: array(mo_tot_num_align,mo_tot_num)
double precision, intent(in) :: value
double precision :: tmp,hij
integer :: i,j
tmp = 0.d0
do i = 1, N_det
do j = 1, N_det
call i_O1_j(array,psi_det(1,1,i),psi_det(1,1,j),N_int,hij)
tmp+= hij * psi_coef(i,1) * psi_coef(j,1)
enddo
enddo
if(dabs(tmp-value).ge.1.d-8)then
print*,'subroutine test_average_value'
print*,'PB WITH AVERAGE VALUE !!!!!!!!!'
print*,'<psi|O|psi> = ',tmp
print*,'value = ',value
stop
endif
end
subroutine test_average_value_alpha_beta(array,value)
implicit none
double precision, intent(in) :: array(mo_tot_num_align,mo_tot_num)
double precision, intent(in) :: value
double precision :: tmp,hij
integer :: i,j
tmp = 0.d0
do i = 1, N_det
do j = 1, N_det
call i_O1_j_alpha_beta(array,psi_det(1,1,i),psi_det(1,1,j),N_int,hij)
tmp+= hij * psi_coef(i,1) * psi_coef(j,1)
enddo
enddo
print*,'tmp = ',tmp
print*,'value = ',value
tmp = 0.d0
do i = 1, N_det
call i_O1_psi_alpha_beta(array,psi_det(1,1,i),psi_det,psi_coef,N_int,N_det,psi_det_size,N_states,hij)
tmp+= hij * psi_coef(i,1)
enddo
print*,'tmp = ',tmp
print*,'value = ',value
if(dabs(tmp-value).ge.1.d-8)then
print*,'subroutine test_average_value'
print*,'PB WITH AVERAGE VALUE !!!!!!!!!'
print*,'<psi|O|psi> = ',tmp
print*,'value = ',value
stop
endif
end
subroutine test_dm(tmp_array)
use bitmasks
implicit none
double precision,intent(out) :: tmp_array(mo_tot_num,mo_tot_num)
double precision :: phase
integer :: i,j,ispin,k
integer :: degree
integer :: exc(0:2,2,2),h1,p1,h2,p2,s1,s2
integer :: occ_det(N_int*bit_kind_size,2)
integer :: tmp
double precision :: accu
tmp_array = 0.d0
do i = 1, N_det
call bitstring_to_list(psi_det(1,1,i), occ_det(1,1), tmp, N_int)
call bitstring_to_list(psi_det(1,2,i), occ_det(1,2), tmp, N_int)
do ispin = 1, 2
do k = 1, elec_num_tab(ispin)
tmp_array(occ_det(k,ispin),occ_det(k,ispin)) += psi_coef(i,1) * psi_coef(i,1)
enddo
enddo
do j = i+1, N_det
call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int)
if (degree.ne.1)cycle
call get_excitation(psi_det(1,1,i),psi_det(1,1,j),exc,degree,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
tmp_array(h1,p1) += phase * psi_coef(i,1) * psi_coef(j,1)
tmp_array(p1,h1) = tmp_array(h1,p1)
enddo
enddo
end

339
src/Properties/slater_rules_mono_electronic.irp.f Normal file
View File

@ -0,0 +1,339 @@
subroutine i_O1_j(array,key_i,key_j,Nint,hij)
use bitmasks
implicit none
BEGIN_DOC
! Returns <i|O1|j> where i and j are determinants
! and O1 is a ONE BODY OPERATOR
! array is the array of the mono electronic operator
! on the MO basis
END_DOC
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2)
double precision, intent(out) :: hij
double precision, intent(in) :: array(mo_tot_num_align,mo_tot_num)
integer :: exc(0:2,2,2)
integer :: degree
integer :: m,p
double precision :: diag_O1_mat_elem, phase
ASSERT (Nint > 0)
ASSERT (Nint == N_int)
ASSERT (sum(popcnt(key_i(:,1))) == elec_alpha_num)
ASSERT (sum(popcnt(key_i(:,2))) == elec_beta_num)
ASSERT (sum(popcnt(key_j(:,1))) == elec_alpha_num)
ASSERT (sum(popcnt(key_j(:,2))) == elec_beta_num)
hij = 0.d0
!DEC$ FORCEINLINE
call get_excitation_degree(key_i,key_j,degree,Nint)
select case (degree)
case (2)
hij = 0.d0
case (1)
call get_mono_excitation(key_i,key_j,exc,phase,Nint)
if (exc(0,1,1) == 1) then
! Mono alpha
m = exc(1,1,1)
p = exc(1,2,1)
else
! Mono beta
m = exc(1,1,2)
p = exc(1,2,2)
endif
hij = phase* array(m,p)
case (0)
hij = diag_O1_mat_elem(array,key_i,Nint)
end select
end
subroutine i_O1_psi(array,key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array)
use bitmasks
implicit none
integer, intent(in) :: Nint, Ndet,Ndet_max,Nstate
double precision, intent(in) :: array(mo_tot_num_align,mo_tot_num)
integer(bit_kind), intent(in) :: keys(Nint,2,Ndet)
integer(bit_kind), intent(in) :: key(Nint,2)
double precision, intent(in) :: coef(Ndet_max,Nstate)
double precision, intent(out) :: i_H_psi_array(Nstate)
integer :: i, ii,j
double precision :: phase
integer :: exc(0:2,2,2)
double precision :: hij
integer :: idx(0:Ndet)
BEGIN_DOC
! <key|O1|psi> for the various Nstates
! and O1 is a ONE BODY OPERATOR
! array is the array of the mono electronic operator
! on the MO basis
END_DOC
ASSERT (Nint > 0)
ASSERT (N_int == Nint)
ASSERT (Nstate > 0)
ASSERT (Ndet > 0)
ASSERT (Ndet_max >= Ndet)
i_H_psi_array = 0.d0
call filter_connected_mono(keys,key,Nint,Ndet,idx)
do ii=1,idx(0)
i = idx(ii)
!DEC$ FORCEINLINE
call i_O1_j(array,keys(1,1,i),key,Nint,hij)
do j = 1, Nstate
i_H_psi_array(j) = i_H_psi_array(j) + coef(i,j)*hij
enddo
enddo
end
double precision function diag_O1_mat_elem(array,det_in,Nint)
use bitmasks
implicit none
BEGIN_DOC
! Computes <i|O1|i>
END_DOC
integer,intent(in) :: Nint
integer(bit_kind),intent(in) :: det_in(Nint,2)
double precision, intent(in) :: array(mo_tot_num_align,mo_tot_num)
integer :: i, ispin,tmp
integer :: occ_det(Nint*bit_kind_size,2)
ASSERT (Nint > 0)
ASSERT (sum(popcnt(det_in(:,1))) == elec_alpha_num)
ASSERT (sum(popcnt(det_in(:,2))) == elec_beta_num)
call bitstring_to_list(det_in(1,1), occ_det(1,1), tmp, Nint)
call bitstring_to_list(det_in(1,2), occ_det(1,2), tmp, Nint)
diag_O1_mat_elem = 0.d0
do ispin = 1, 2
do i = 1, elec_num_tab(ispin)
diag_O1_mat_elem += array(occ_det(i,ispin),occ_det(i,ispin))
enddo
enddo
end
subroutine i_O1_psi_alpha_beta(array,key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array)
use bitmasks
implicit none
integer, intent(in) :: Nint, Ndet,Ndet_max,Nstate
double precision, intent(in) :: array(mo_tot_num_align,mo_tot_num)
integer(bit_kind), intent(in) :: keys(Nint,2,Ndet)
integer(bit_kind), intent(in) :: key(Nint,2)
double precision, intent(in) :: coef(Ndet_max,Nstate)
double precision, intent(out) :: i_H_psi_array(Nstate)
integer :: i, ii,j
double precision :: phase
integer :: exc(0:2,2,2)
double precision :: hij
integer :: idx(0:Ndet)
BEGIN_DOC
! <key|O1(alpha) - O1(beta)|psi> for the various Nstates
! and O1 is a ONE BODY OPERATOR
! array is the array of the mono electronic operator
! on the MO basis
END_DOC
ASSERT (Nint > 0)
ASSERT (N_int == Nint)
ASSERT (Nstate > 0)
ASSERT (Ndet > 0)
ASSERT (Ndet_max >= Ndet)
i_H_psi_array = 0.d0
call filter_connected_mono(keys,key,Nint,Ndet,idx)
do ii=1,idx(0)
i = idx(ii)
!DEC$ FORCEINLINE
call i_O1_j_alpha_beta(array,keys(1,1,i),key,Nint,hij)
do j = 1, Nstate
i_H_psi_array(j) = i_H_psi_array(j) + coef(i,j)*hij
enddo
enddo
end
subroutine i_O1_j_alpha_beta(array,key_i,key_j,Nint,hij)
use bitmasks
implicit none
BEGIN_DOC
! Returns <i|O1(alpha) - O1(beta)|j> where i and j are determinants
! and O1 is a ONE BODY OPERATOR
! array is the array of the mono electronic operator
! on the MO basis
END_DOC
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2)
double precision, intent(out) :: hij
double precision, intent(in) :: array(mo_tot_num_align,mo_tot_num)
integer :: exc(0:2,2,2)
integer :: degree
integer :: m,p
double precision :: diag_O1_mat_elem_alpha_beta, phase
ASSERT (Nint > 0)
ASSERT (Nint == N_int)
ASSERT (sum(popcnt(key_i(:,1))) == elec_alpha_num)
ASSERT (sum(popcnt(key_i(:,2))) == elec_beta_num)
ASSERT (sum(popcnt(key_j(:,1))) == elec_alpha_num)
ASSERT (sum(popcnt(key_j(:,2))) == elec_beta_num)
hij = 0.d0
!DEC$ FORCEINLINE
call get_excitation_degree(key_i,key_j,degree,Nint)
select case (degree)
case (2)
hij = 0.d0
case (1)
call get_mono_excitation(key_i,key_j,exc,phase,Nint)
if (exc(0,1,1) == 1) then
! Mono alpha
m = exc(1,1,1)
p = exc(1,2,1)
hij = phase* array(m,p)
else
! Mono beta
m = exc(1,1,2)
p = exc(1,2,2)
hij = -phase* array(m,p)
endif
case (0)
hij = diag_O1_mat_elem_alpha_beta(array,key_i,Nint)
end select
end
double precision function diag_O1_mat_elem_alpha_beta(array,det_in,Nint)
use bitmasks
implicit none
BEGIN_DOC
! Computes <i|O1(alpha) -O1(beta)|i>
END_DOC
integer,intent(in) :: Nint
integer(bit_kind),intent(in) :: det_in(Nint,2)
double precision, intent(in) :: array(mo_tot_num_align,mo_tot_num)
integer :: i, ispin,tmp
integer :: occ_det(Nint*bit_kind_size,2)
ASSERT (Nint > 0)
ASSERT (sum(popcnt(det_in(:,1))) == elec_alpha_num)
ASSERT (sum(popcnt(det_in(:,2))) == elec_beta_num)
call bitstring_to_list(det_in(1,1), occ_det(1,1), tmp, Nint)
call bitstring_to_list(det_in(1,2), occ_det(1,2), tmp, Nint)
diag_O1_mat_elem_alpha_beta = 0.d0
ispin = 1
do i = 1, elec_num_tab(ispin)
diag_O1_mat_elem_alpha_beta += array(occ_det(i,ispin),occ_det(i,ispin))
enddo
ispin = 2
do i = 1, elec_num_tab(ispin)
diag_O1_mat_elem_alpha_beta -= array(occ_det(i,ispin),occ_det(i,ispin))
enddo
end
subroutine filter_connected_mono(key1,key2,Nint,sze,idx)
use bitmasks
implicit none
BEGIN_DOC
! Filters out the determinants that are not connected through PURE
!
! MONO EXCITATIONS OPERATORS (a^{\dagger}j a_i)
!
! returns the array idx which contains the index of the
!
! determinants in the array key1 that interact
!
! via some PURE MONO EXCITATIONS OPERATORS
!
! idx(0) is the number of determinants that interact with key1
END_DOC
integer, intent(in) :: Nint, sze
integer(bit_kind), intent(in) :: key1(Nint,2,sze)
integer(bit_kind), intent(in) :: key2(Nint,2)
integer, intent(out) :: idx(0:sze)
integer :: i,j,l
integer :: degree_x2
ASSERT (Nint > 0)
ASSERT (sze >= 0)
l=1
if (Nint==1) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = popcnt( xor( key1(1,1,i), key2(1,1))) &
+ popcnt( xor( key1(1,2,i), key2(1,2)))
if (degree_x2 > 3) then
cycle
else
idx(l) = i
l = l+1
endif
enddo
else if (Nint==2) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + &
popcnt(xor( key1(2,1,i), key2(2,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2))) + &
popcnt(xor( key1(2,2,i), key2(2,2)))
if (degree_x2 > 3) then
cycle
else
idx(l) = i
l = l+1
endif
enddo
else if (Nint==3) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2))) + &
popcnt(xor( key1(2,1,i), key2(2,1))) + &
popcnt(xor( key1(2,2,i), key2(2,2))) + &
popcnt(xor( key1(3,1,i), key2(3,1))) + &
popcnt(xor( key1(3,2,i), key2(3,2)))
if (degree_x2 > 3) then
cycle
else
idx(l) = i
l = l+1
endif
enddo
else
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = 0
!DEC$ LOOP COUNT MIN(4)
do j=1,Nint
degree_x2 = degree_x2+ popcnt(xor( key1(j,1,i), key2(j,1))) +&
popcnt(xor( key1(j,2,i), key2(j,2)))
if (degree_x2 > 3) then
exit
endif
enddo
if (degree_x2 <= 3) then
idx(l) = i
l = l+1
endif
enddo
endif
idx(0) = l-1
end

12
src/Selectors_full/selectors.irp.f
View File

@ -60,6 +60,18 @@ END_PROVIDER
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, psi_selectors_diag_h_mat, (psi_selectors_size) ]
implicit none
BEGIN_DOC
! Diagonal elements of the H matrix for each selectors
END_DOC
integer :: i
double precision :: diag_H_mat_elem
do i = 1, N_det_selectors
psi_selectors_diag_h_mat(i) = diag_H_mat_elem(psi_selectors(1,1,i),N_int)
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), psi_selectors_ab, (N_int,2,psi_selectors_size) ]
&BEGIN_PROVIDER [ double precision, psi_selectors_coef_ab, (psi_selectors_size,N_states) ]

0
src/Selectors_no_sorted/ASSUMPTIONS.rst Normal file
View File

8
src/Selectors_no_sorted/Makefile Normal file
View File

@ -0,0 +1,8 @@
default: all
# Define here all new external source files and objects.Don't forget to prefix the
# object files with IRPF90_temp/
SRC=
OBJ=
include $(QPACKAGE_ROOT)/src/Makefile.common

1
src/Selectors_no_sorted/NEEDED_MODULES Normal file
View File

@ -0,0 +1 @@
AOs BiInts Bitmask Dets Electrons Ezfio_files MonoInts MOs Nuclei Output Utils

4
src/Selectors_no_sorted/README.rst Normal file
View File

@ -0,0 +1,4 @@
==========================
Selectors_no_sorted Module
==========================

79
src/Selectors_no_sorted/e_corr_selectors.irp.f Normal file
View File

@ -0,0 +1,79 @@
use bitmasks
BEGIN_PROVIDER [integer, exc_degree_per_selectors, (N_det_selectors)]
&BEGIN_PROVIDER [integer, double_index_selectors, (N_det_selectors)]
&BEGIN_PROVIDER [integer, n_double_selectors]
implicit none
BEGIN_DOC
! degree of excitation respect to Hartree Fock for the wave function
!
! for the all the selectors determinants
!
! double_index_selectors = list of the index of the double excitations
!
! n_double_selectors = number of double excitations in the selectors determinants
END_DOC
integer :: i,degree
n_double_selectors = 0
do i = 1, N_det_selectors
call get_excitation_degree(psi_selectors(1,1,i),ref_bitmask,degree,N_int)
exc_degree_per_selectors(i) = degree
if(degree==2)then
n_double_selectors += 1
double_index_selectors(n_double_selectors) =i
endif
enddo
END_PROVIDER
BEGIN_PROVIDER[double precision, coef_hf_selector]
&BEGIN_PROVIDER[double precision, inv_selectors_coef_hf]
&BEGIN_PROVIDER[double precision, inv_selectors_coef_hf_squared]
&BEGIN_PROVIDER[double precision, E_corr_per_selectors, (N_det_selectors)]
&BEGIN_PROVIDER[double precision, i_H_HF_per_selectors, (N_det_selectors)]
&BEGIN_PROVIDER[double precision, Delta_E_per_selector, (N_det_selectors)]
&BEGIN_PROVIDER[double precision, E_corr_double_only ]
&BEGIN_PROVIDER[double precision, E_corr_second_order ]
implicit none
BEGIN_DOC
! energy of correlation per determinant respect to the Hartree Fock determinant
!
! for the all the double excitations in the selectors determinants
!
! E_corr_per_selectors(i) = <D_i|H|HF> * c(D_i)/c(HF) if |D_i> is a double excitation
!
! E_corr_per_selectors(i) = -1000.d0 if it is not a double excitation
!
! coef_hf_selector = coefficient of the Hartree Fock determinant in the selectors determinants
END_DOC
PROVIDE ref_bitmask_energy psi_selectors ref_bitmask N_int psi_selectors
integer :: i,degree
double precision :: hij,diag_H_mat_elem
E_corr_double_only = 0.d0
E_corr_second_order = 0.d0
do i = 1, N_det_selectors
if(exc_degree_per_selectors(i)==2)then
call i_H_j(ref_bitmask,psi_selectors(1,1,i),N_int,hij)
i_H_HF_per_selectors(i) = hij
E_corr_per_selectors(i) = psi_selectors_coef(i,1) * hij
E_corr_double_only += E_corr_per_selectors(i)
E_corr_second_order += hij * hij /(ref_bitmask_energy - diag_H_mat_elem(psi_selectors(1,1,i),N_int))
elseif(exc_degree_per_selectors(i) == 0)then
coef_hf_selector = psi_selectors_coef(i,1)
E_corr_per_selectors(i) = -1000.d0
Delta_E_per_selector(i) = 0.d0
else
E_corr_per_selectors(i) = -1000.d0
endif
enddo
if (dabs(coef_hf_selector) > 1.d-8) then
inv_selectors_coef_hf = 1.d0/coef_hf_selector
inv_selectors_coef_hf_squared = inv_selectors_coef_hf * inv_selectors_coef_hf
else
inv_selectors_coef_hf = 0.d0
inv_selectors_coef_hf_squared = 0.d0
endif
do i = 1,n_double_selectors
E_corr_per_selectors(double_index_selectors(i)) *=inv_selectors_coef_hf
enddo
E_corr_double_only = E_corr_double_only * inv_selectors_coef_hf
END_PROVIDER

156
src/Selectors_no_sorted/selectors.irp.f Normal file
View File

@ -0,0 +1,156 @@
use bitmasks
BEGIN_PROVIDER [ integer, psi_selectors_size ]
implicit none
psi_selectors_size = psi_det_size
END_PROVIDER
BEGIN_PROVIDER [ integer, N_det_selectors]
implicit none
BEGIN_DOC
! For Single reference wave functions, the number of selectors is 1 : the
! Hartree-Fock determinant
END_DOC
integer :: i
double precision :: norm
call write_time(output_dets)
norm = 0.d0
N_det_selectors = N_det
N_det_selectors = max(N_det_selectors,1)
call write_int(output_dets,N_det_selectors,'Number of selectors')
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), psi_selectors, (N_int,2,psi_selectors_size) ]
&BEGIN_PROVIDER [ double precision, psi_selectors_coef, (psi_selectors_size,N_states) ]
implicit none
BEGIN_DOC
! Determinants on which we apply <i|H|psi> for perturbation.
END_DOC
integer :: i,k
do i=1,N_det_selectors
do k=1,N_int
psi_selectors(k,1,i) = psi_det(k,1,i)
psi_selectors(k,2,i) = psi_det(k,2,i)
enddo
enddo
do k=1,N_states
do i=1,N_det_selectors
psi_selectors_coef(i,k) = psi_coef(i,k)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, psi_selectors_diag_h_mat, (psi_selectors_size) ]
implicit none
BEGIN_DOC
! Diagonal elements of the H matrix for each selectors
END_DOC
integer :: i
double precision :: diag_H_mat_elem
do i = 1, N_det_selectors
psi_selectors_diag_h_mat(i) = diag_H_mat_elem(psi_selectors(1,1,i),N_int)
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), psi_selectors_ab, (N_int,2,psi_selectors_size) ]
&BEGIN_PROVIDER [ double precision, psi_selectors_coef_ab, (psi_selectors_size,N_states) ]
&BEGIN_PROVIDER [ integer, psi_selectors_next_ab, (2,psi_selectors_size) ]
implicit none
BEGIN_DOC
! Determinants on which we apply <i|H|j>.
! They are sorted by the 3 highest electrons in the alpha part,
! then by the 3 highest electrons in the beta part to accelerate
! the research of connected determinants.
END_DOC
integer :: i,j,k
integer, allocatable :: iorder(:)
integer*8, allocatable :: bit_tmp(:)
integer*8, external :: det_search_key
allocate ( iorder(N_det_selectors), bit_tmp(N_det_selectors) )
! Sort alpha dets
! ---------------
integer(bit_kind) :: det_tmp(N_int)
do i=1,N_det_selectors
iorder(i) = i
call int_of_3_highest_electrons(psi_selectors(1,1,i),bit_tmp(i),N_int)
enddo
call i8sort(bit_tmp,iorder,N_det_selectors)
!DIR$ IVDEP
do i=1,N_det_selectors
do j=1,N_int
psi_selectors_ab(j,1,i) = psi_selectors(j,1,iorder(i))
psi_selectors_ab(j,2,i) = psi_selectors(j,2,iorder(i))
enddo
do k=1,N_states
psi_coef_sorted_ab(i,k) = psi_selectors_coef(iorder(i),k)
enddo
enddo
! Find next alpha
! ---------------
integer :: next
next = N_det_selectors+1
psi_selectors_next_ab(1,N_det_selectors) = next
do i=N_det_selectors-1,1,-1
if (bit_tmp(i) /= bit_tmp(i+1)) then
next = i+1
endif
psi_selectors_next_ab(1,i) = next
enddo
! Sort beta dets
! --------------
integer :: istart, iend
integer(bit_kind), allocatable :: psi_selectors_ab_temp (:,:)
allocate ( psi_selectors_ab_temp (N_int,N_det_selectors) )
do i=1,N_det_selectors
do j=1,N_int
psi_selectors_ab_temp(j,i) = psi_selectors_ab(j,2,i)
enddo
iorder(i) = i
call int_of_3_highest_electrons(psi_selectors_ab_temp(1,i),bit_tmp(i),N_int)
enddo
istart=1
do while ( istart<N_det_selectors )
iend = psi_selectors_next_ab(1,istart)
call i8sort(bit_tmp(istart),iorder(istart),iend-istart)
!DIR$ IVDEP
do i=istart,iend-1
do j=1,N_int
psi_selectors_ab(j,2,i) = psi_selectors_ab_temp(j,iorder(i))
enddo
do k=1,N_states
psi_coef_sorted_ab(i,k) = psi_coef(iorder(i),k)
enddo
enddo
next = iend
psi_selectors_next_ab(2,iend-1) = next
do i=iend-2,1,-1
if (bit_tmp(i) /= bit_tmp(i+1)) then
next = i+1
endif
psi_selectors_next_ab(2,i) = next
enddo
istart = iend
enddo
deallocate(iorder, bit_tmp, psi_selectors_ab_temp)
END_PROVIDER