LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-25 22:52:15 +02:00
Code Issues Releases Wiki Activity

CISD_SC2 ok, CISD_SC2_selected ok

Browse Source
This commit is contained in:
Manu 2014-05-30 18:07:04 +02:00
parent cbc70ec1db
commit 868ef807bd
30 changed files with 532 additions and 302 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
scripts/generate_h_apply.py
View File

@ -123,7 +123,8 @@ class H_apply(object):
sum_norm_pert(k) = 0.d0
sum_H_pert_diag(k) = 0.d0
enddo
"""
"""
self.data["deinit_thread"] = """
!$OMP CRITICAL
do k=1,N_st
@ -136,7 +137,7 @@ class H_apply(object):
"""
self.data["size_max"] = "256"
self.data["initialization"] = """
PROVIDE CI_electronic_energy psi_selectors_coef psi_selectors
PROVIDE CI_electronic_energy psi_selectors_coef psi_selectors E_corr_per_selectors
"""
self.data["keys_work"] = """
call perturb_buffer_%s(i_generator,keys_out,key_idx,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert, &

21
src/CISD/README.rst
View File

@ -35,26 +35,7 @@ Documentation
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
`cisd_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/CISD/SC2.irp.f#L1>`_
CISD+SC2 method :: take off all the disconnected terms of a CISD (selected or not)
.br
dets_in : bitmasks corresponding to determinants
.br
u_in : guess coefficients on the various states. Overwritten
on exit
.br
dim_in : leftmost dimension of u_in
.br
sze : Number of determinants
.br
N_st : Number of eigenstates
.br
Initial guess vectors are not necessarily orthonormal
`repeat_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/CISD/SC2.irp.f#L143>`_
Undocumented
`cisd <http://github.com/LCPQ/quantum_package/tree/master/src/CISD/cisd_sc2.irp.f#L1>`_
`cisd <http://github.com/LCPQ/quantum_package/tree/master/src/CISD/cisd_lapack.irp.f#L1>`_
Undocumented

7
src/CISD/cisd_lapack.irp.f
View File

@ -2,9 +2,8 @@ program cisd
implicit none
integer :: i
N_states = 3
diag_algorithm = "Lapack"
touch N_states diag_algorithm
touch diag_algorithm
print *, 'HF = ', HF_energy
print *, 'N_states = ', N_states
call H_apply_cisd
@ -14,8 +13,4 @@ program cisd
print *, 'E_corr = ',CI_electronic_energy(i) - ref_bitmask_energy
enddo
! call CISD_SC2(psi_det,psi_coef,eigvalues,size(psi_coef,1),N_det,N_states,N_int)
! do i = 1, N_states
! print*,'eigvalues(i) = ',eigvalues(i)
! enddo
end

23
src/CISD/cisd_sc2.irp.f
View File

@ -1,23 +0,0 @@
program cisd
implicit none
integer :: i, j
print *, 'HF = ', HF_energy
print *, 'N_states = ', N_states
call H_apply_cisd
print *, 'N_det = ', N_det
do i = 1,N_states
print *, 'energy = ',CI_energy(i)
print *, 'E_corr = ',CI_electronic_energy(i) - ref_bitmask_energy
do j=1,N_det
psi_coef(j,i) = CI_eigenvectors(j,i)
enddo
enddo
SOFT_TOUCH CI_electronic_energy CI_eigenvectors
call CISD_SC2(psi_det,psi_coef,CI_electronic_energy,size(psi_coef,1),N_det,N_states,N_int)
TOUCH CI_electronic_energy CI_eigenvectors
do i = 1, N_states
print*,'eigvalues(i) = ',CI_energy(i)
enddo
end

34
src/CISD_SC2/README.rst Normal file
View File

@ -0,0 +1,34 @@
===============
CISD_SC2 Module
===============
Documentation
=============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
Needed Modules
==============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
* `AOs <http://github.com/LCPQ/quantum_package/tree/master/src/AOs>`_
* `BiInts <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts>`_
* `Bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask>`_
* `CISD <http://github.com/LCPQ/quantum_package/tree/master/src/CISD>`_
* `Dets <http://github.com/LCPQ/quantum_package/tree/master/src/Dets>`_
* `Electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Electrons>`_
* `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
* `Hartree_Fock <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock>`_
* `MonoInts <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts>`_
* `MOs <http://github.com/LCPQ/quantum_package/tree/master/src/MOs>`_
* `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_
* `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
* `SingleRefMethod <http://github.com/LCPQ/quantum_package/tree/master/src/SingleRefMethod>`_
* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
* `Selectors_full <http://github.com/LCPQ/quantum_package/tree/master/src/Selectors_full>`_

63
src/CISD/SC2.irp.f → src/CISD_SC2/SC2.irp.f
View File

@ -1,4 +1,4 @@
subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint)
subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint,iunit)
use bitmasks
implicit none
BEGIN_DOC
@ -17,7 +17,7 @@ subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint)
!
! Initial guess vectors are not necessarily orthonormal
END_DOC
integer, intent(in) :: dim_in, sze, N_st, Nint
integer, intent(in) :: dim_in, sze, N_st, Nint,iunit
integer(bit_kind), intent(in) :: dets_in(Nint,2,sze)
double precision, intent(inout) :: u_in(dim_in,N_st)
double precision, intent(out) :: energies(N_st)
@ -38,11 +38,22 @@ subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint)
double precision :: e_corr_double_before,accu,cpu_2,cpu_1
integer :: degree_exc(sze)
integer :: i_ok
double precision, allocatable :: eigenvectors(:,:), eigenvalues(:),H_matrix_tmp(:,:)
if(sze<500)then
allocate (eigenvectors(size(H_matrix_all_dets,1),N_det))
allocate (H_matrix_tmp(size(H_matrix_all_dets,1),N_det))
allocate (eigenvalues(N_det))
do i = 1, sze
do j = 1, sze
H_matrix_tmp(i,j) = H_matrix_all_dets(i,j)
enddo
enddo
endif
call write_time(output_CISD)
write(output_CISD,'(A)') ''
write(output_CISD,'(A)') 'CISD SC2'
write(output_CISD,'(A)') '========'
call write_time(output_CISD_SC2)
write(output_CISD_SC2,'(A)') ''
write(output_CISD_SC2,'(A)') 'CISD SC2'
write(output_CISD_SC2,'(A)') '========'
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP SHARED(sze,N_st, &
!$OMP H_jj_ref,Nint,dets_in,u_in) &
@ -108,25 +119,40 @@ subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint)
H_jj_dressed(i) += accu
enddo
call davidson_diag_hjj(dets_in,u_in,H_jj_dressed,energies,dim_in,sze,N_st,Nint,output_CISD)
if(sze>500)then
call davidson_diag_hjj(dets_in,u_in,H_jj_dressed,energies,dim_in,sze,N_st,Nint,output_CISD_SC2)
else
do i = 1,sze
H_matrix_tmp(i,i) = H_jj_dressed(i)
enddo
call lapack_diag(eigenvalues,eigenvectors, &
H_matrix_tmp,size(H_matrix_all_dets,1),N_det)
do j=1,min(N_states,sze)
do i=1,sze
u_in(i,j) = eigenvectors(i,j)
enddo
energies(j) = eigenvalues(j)
enddo
endif
e_corr_double = 0.d0
inv_c0 = 1.d0/u_in(index_hf,1)
do i = 1, N_double
e_corr_array(i) = u_in(index_double(i),1)*inv_c0 * hij_double(i)
e_corr_double += e_corr_array(i)
enddo
write(output_CISD,'(A,I3)') 'SC2 Iteration ', iter
write(output_CISD,'(A)') '------------------'
write(output_CISD,'(A)') ''
write(output_CISD,'(A)') '===== ================'
write(output_CISD,'(A)') 'State Energy '
write(output_CISD,'(A)') '===== ================'
write(output_CISD_SC2,'(A,I3)') 'SC2 Iteration ', iter
write(output_CISD_SC2,'(A)') '------------------'
write(output_CISD_SC2,'(A)') ''
write(output_CISD_SC2,'(A)') '===== ================'
write(output_CISD_SC2,'(A)') 'State Energy '
write(output_CISD_SC2,'(A)') '===== ================'
do i=1,N_st
write(output_CISD,'(I5,X,F16.10)') i, energies(i)+nuclear_repulsion
write(output_CISD_SC2,'(I5,X,F16.10)') i, energies(i)+nuclear_repulsion
enddo
write(output_CISD,'(A)') '===== ================'
write(output_CISD,'(A)') ''
call write_double(output_CISD,(e_corr_double - e_corr_double_before),&
write(output_CISD_SC2,'(A)') '===== ================'
write(output_CISD_SC2,'(A)') ''
call write_double(output_CISD_SC2,(e_corr_double - e_corr_double_before),&
'Delta(E_corr)')
converged = dabs(e_corr_double - e_corr_double_before) < 1.d-10
if (converged) then
@ -136,7 +162,7 @@ subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint)
enddo
call write_time(output_CISD)
call write_time(output_CISD_SC2)
end
@ -187,3 +213,4 @@ subroutine repeat_excitation(key_in,key_1,key_2,i_ok,Nint)
return
endif
end

14
src/CISD_SC2/cisd_SC2.irp.f Normal file
View File

@ -0,0 +1,14 @@
program cisd
implicit none
integer :: i
print *, ' HF = ', HF_energy
print *, 'N_states = ', N_states
call H_apply_cisd
print *, 'N_det = ', N_det
do i = 1,N_states
print *, 'energy = ',CI_SC2_energy(i)
print *, 'E_corr = ',CI_SC2_electronic_energy(i) - ref_bitmask_energy
enddo
end

69
src/CISD_SC2/diagonalize_CI_SC2.irp.f Normal file
View File

@ -0,0 +1,69 @@
BEGIN_PROVIDER [ character*(64), diag_algorithm ]
implicit none
BEGIN_DOC
! Diagonalization algorithm (Davidson or Lapack)
END_DOC
if (N_det > 500) then
diag_algorithm = "Davidson"
else
diag_algorithm = "Lapack"
endif
if (N_det < N_states) then
diag_algorithm = "Lapack"
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, CI_SC2_energy, (N_states) ]
implicit none
BEGIN_DOC
! N_states lowest eigenvalues of the CI matrix
END_DOC
integer :: j
character*(8) :: st
call write_time(output_CISD_SC2)
do j=1,N_states
CI_SC2_energy(j) = CI_SC2_electronic_energy(j) + nuclear_repulsion
write(st,'(I4)') j
call write_double(output_CISD_SC2,CI_SC2_energy(j),'Energy of state '//trim(st))
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, CI_SC2_electronic_energy, (N_states) ]
&BEGIN_PROVIDER [ double precision, CI_SC2_eigenvectors, (N_det,N_states) ]
implicit none
BEGIN_DOC
! Eigenvectors/values of the CI matrix
END_DOC
integer :: i,j
do j=1,N_states
do i=1,N_det
CI_SC2_eigenvectors(i,j) = CI_eigenvectors(i,j)
enddo
CI_SC2_electronic_energy(j) = CI_electronic_energy(j)
enddo
print*,'E(CISD) = ',CI_electronic_energy(1)
call CISD_SC2(psi_det,CI_SC2_eigenvectors,CI_SC2_electronic_energy, &
size(CI_SC2_eigenvectors,1),N_det,N_states,N_int,output_CISD_SC2)
END_PROVIDER
subroutine diagonalize_CI_SC2
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
do i=1,N_det
psi_coef(i,j) = CI_SC2_eigenvectors(i,j)
enddo
enddo
SOFT_TOUCH psi_coef psi_det CI_SC2_electronic_energy CI_SC2_energy CI_SC2_eigenvectors
end

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

8
src/CISD_SC2_selected/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/CISD_SC2_selected/NEEDED_MODULES Normal file
View File

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

37
src/CISD_SC2_selected/README.rst Normal file
View File

@ -0,0 +1,37 @@
========================
CISD_SC2_selected Module
========================
Documentation
=============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
Needed Modules
==============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
* `AOs <http://github.com/LCPQ/quantum_package/tree/master/src/AOs>`_
* `BiInts <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts>`_
* `Bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask>`_
* `CISD <http://github.com/LCPQ/quantum_package/tree/master/src/CISD>`_
* `CISD_SC2 <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_SC2>`_
* `CISD_selected <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_selected>`_
* `Dets <http://github.com/LCPQ/quantum_package/tree/master/src/Dets>`_
* `Electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Electrons>`_
* `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
* `Hartree_Fock <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock>`_
* `MonoInts <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts>`_
* `MOs <http://github.com/LCPQ/quantum_package/tree/master/src/MOs>`_
* `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_
* `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
* `Perturbation <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation>`_
* `Selectors_full <http://github.com/LCPQ/quantum_package/tree/master/src/Selectors_full>`_
* `SingleRefMethod <http://github.com/LCPQ/quantum_package/tree/master/src/SingleRefMethod>`_
* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_

32
src/CISD_SC2_selected/cisd_sc2_selection.irp.f Normal file
View File

@ -0,0 +1,32 @@
program cisd_sc2_selected
implicit none
integer :: i,k
double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:),E_old(:)
integer :: N_st, iter
character*(64) :: perturbation
N_st = N_states
allocate (pt2(N_st), norm_pert(N_st), H_pert_diag(N_st),E_old(N_st))
pt2 = 1.d0
perturbation = "epstein_nesbet_sc2_projected"
E_old(1) = HF_energy
do while (maxval(abs(pt2(1:N_st))) > 1.d-9)
print*,'----'
print*,''
call H_apply_cisd_selection(perturbation,pt2, norm_pert, H_pert_diag, N_st)
call diagonalize_CI_SC2
print *, 'N_det = ', N_det
print *, 'PT2(SC2) = ', pt2
print *, 'E(SC2) = ', CI_SC2_energy(1)
print *, 'E_before(SC2)+PT2(SC2) = ', (E_old(1)+pt2(1))
print *, 'E(SC2)+PT2(projctd)SC2 = ', (E_old(1)+H_pert_diag(1))
! print *, 'E corr = ', (E_old(1)) - HF_energy
E_old(1) = CI_SC2_energy(1)
if (abort_all) then
exit
endif
enddo
deallocate(pt2,norm_pert,H_pert_diag)
end

2
src/CISD_selected/NEEDED_MODULES
View File

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

1
src/CISD_selected/README.rst
View File

@ -26,6 +26,7 @@ Needed Modules
* `BiInts <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts>`_
* `Bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask>`_
* `CISD <http://github.com/LCPQ/quantum_package/tree/master/src/CISD>`_
* `CISD_SC2 <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_SC2>`_
* `Dets <http://github.com/LCPQ/quantum_package/tree/master/src/Dets>`_
* `Electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Electrons>`_
* `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_

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

@ -3,22 +3,27 @@ program cisd
integer :: i,k
double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:)
double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:),E_old(:)
integer :: N_st, iter
character*(64) :: perturbation
N_st = N_states
allocate (pt2(N_st), norm_pert(N_st), H_pert_diag(N_st))
allocate (pt2(N_st), norm_pert(N_st), H_pert_diag(N_st),E_old(N_st))
pt2 = 1.d0
perturbation = "epstein_nesbet"
E_old(1) = HF_energy
do while (maxval(abs(pt2(1:N_st))) > 1.d-6)
print*,'----'
print*,''
call H_apply_cisd_selection(perturbation,pt2, norm_pert, H_pert_diag, N_st)
call diagonalize_CI
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
print *, 'PT2 = ', pt2
print *, 'E_old = ', E_old(1)
print *, 'E = ', CI_energy
print *, 'E+PT2 = ', CI_energy+pt2
print *, 'E+PT2 = ', E_old(1)+pt2(1)
! print *, 'E+PT2_new= ', (E_old(1)+1.d0*pt2(1)+H_pert_diag(1))/(1.d0 +norm_pert(1))
E_old(1) = CI_energy(1)
if (abort_all) then
exit
endif

4
src/Dets/Makefile
View File

@ -2,11 +2,7 @@ default: all
# Define here all new external source files and objects.Don't forget to prefix the
# object files with IRPF90_temp/
<<<<<<< HEAD
SRC=
=======
SRC=H_apply_template.f
>>>>>>> 186b3a9e572d016a28ee0cbc3caf0be88a290931
OBJ=
include $(QPACKAGE_ROOT)/src/Makefile.common

64
src/Dets/README.rst
View File

@ -51,9 +51,10 @@ Documentation
.. NEEDED_MODULES file.
`copy_h_apply_buffer_to_wf <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L113>`_
Undocumented
Copies the H_apply buffer to psi_coef. You need to touch psi_det, psi_coef and N_det
after calling this function.
`fill_h_apply_buffer_no_selection <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L197>`_
`fill_h_apply_buffer_no_selection <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L199>`_
Fill the H_apply buffer with determiants for CISD
`h_apply_buffer_allocated <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L14>`_
@ -69,18 +70,18 @@ Documentation
`connected_to_ref <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/connected_to_ref.irp.f#L1>`_
Undocumented
`det_is_not_or_may_be_in_ref <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/connected_to_ref.irp.f#L196>`_
`det_is_not_or_may_be_in_ref <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/connected_to_ref.irp.f#L188>`_
If true, det is not in ref
If false, det may be in ref
`key_pattern_not_in_ref <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/connected_to_ref.irp.f#L229>`_
`key_pattern_not_in_ref <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/connected_to_ref.irp.f#L222>`_
Min and max values of the integers of the keys of the reference
`davidson_converged <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L375>`_
`davidson_converged <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L381>`_
True if the Davidson algorithm is converged
`davidson_criterion <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L365>`_
Can be : [ energy | residual | both ]
`davidson_criterion <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L371>`_
Can be : [ energy | residual | both | wall_time | cpu_time | iterations ]
`davidson_diag <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L18>`_
Davidson diagonalization.
@ -126,24 +127,41 @@ Documentation
`davidson_sze_max <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L9>`_
Max number of Davidson sizes
`davidson_threshold <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L366>`_
Can be : [ energy | residual | both ]
`davidson_threshold <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L372>`_
Can be : [ energy | residual | both | wall_time | cpu_time | iterations ]
`n_det <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L11>`_
`n_det <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L20>`_
Number of determinants in the wave function
`n_det_reference <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L75>`_
Number of determinants in the reference wave function
`n_states <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L3>`_
Number of states to consider
`psi_coef <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L28>`_
The wave function. Initialized with Hartree-Fock
`psi_average_norm_contrib <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L84>`_
Contribution of determinants to the state-averaged density
`psi_det <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L27>`_
The wave function. Initialized with Hartree-Fock
`psi_average_norm_contrib_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L105>`_
Wave function sorted by determinants (state-averaged)
`psi_det_size <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L19>`_
`psi_coef <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L47>`_
The wave function. Initialized with Hartree-Fock if the EZFIO file
is empty
`psi_coef_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L104>`_
Wave function sorted by determinants (state-averaged)
`psi_det <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L46>`_
The wave function. Initialized with Hartree-Fock if the EZFIO file
is empty
`psi_det_size <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L38>`_
Size of the psi_det/psi_coef arrays
`psi_det_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L103>`_
Wave function sorted by determinants (state-averaged)
`double_exc_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants_bitmasks.irp.f#L40>`_
double_exc_bitmask(:,1,i) is the bitmask for holes of excitation 1
double_exc_bitmask(:,2,i) is the bitmask for particles of excitation 1
@ -162,6 +180,22 @@ Documentation
single_exc_bitmask(:,2,i) is the bitmask for particles
for a given couple of hole/particle excitations i.
`ci_eigenvectors <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/diagonalize_CI.irp.f#L36>`_
Eigenvectors/values of the CI matrix
`ci_electronic_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/diagonalize_CI.irp.f#L35>`_
Eigenvectors/values of the CI matrix
`ci_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/diagonalize_CI.irp.f#L18>`_
N_states lowest eigenvalues of the CI matrix
`diag_algorithm <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/diagonalize_CI.irp.f#L1>`_
Diagonalization algorithm (Davidson or Lapack)
`diagonalize_ci <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/diagonalize_CI.irp.f#L73>`_
Replace the coefficients of the CI states by the coefficients of the
eigenstates of the CI matrix
`filter_connected <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/filter_connected.irp.f#L2>`_
Filters out the determinants that are not connected by H
.br

21
src/Dets/filter_connected.irp.f
View File

@ -221,6 +221,7 @@ subroutine filter_connected_i_H_psi0_SC2(key1,key2,Nint,sze,idx,idx_repeat)
l_repeat=1
call get_excitation_degree(ref_bitmask,key2,degree,Nint)
integer :: degree
ASSERT (degree .ne. 0)
if(degree == 2)then
if (Nint==1) then
@ -235,7 +236,7 @@ subroutine filter_connected_i_H_psi0_SC2(key1,key2,Nint,sze,idx,idx_repeat)
idx(l) = i
l = l+1
endif
elseif(degree>6)then
elseif(degree_x2>6)then
idx_repeat(l_repeat) = i
l_repeat = l_repeat + 1
endif
@ -254,7 +255,7 @@ subroutine filter_connected_i_H_psi0_SC2(key1,key2,Nint,sze,idx,idx_repeat)
idx(l) = i
l = l+1
endif
elseif(degree>6)then
elseif(degree_x2>6)then
idx_repeat(l_repeat) = i
l_repeat = l_repeat + 1
endif
@ -275,7 +276,7 @@ subroutine filter_connected_i_H_psi0_SC2(key1,key2,Nint,sze,idx,idx_repeat)
idx(l) = i
l = l+1
endif
elseif(degree>6)then
elseif(degree_x2>6)then
idx_repeat(l_repeat) = i
l_repeat = l_repeat + 1
endif
@ -299,7 +300,7 @@ subroutine filter_connected_i_H_psi0_SC2(key1,key2,Nint,sze,idx,idx_repeat)
idx(l) = i
l = l+1
endif
elseif(degree>6)then
elseif(degree_x2>6)then
idx_repeat(l_repeat) = i
l_repeat = l_repeat + 1
endif
@ -392,11 +393,13 @@ subroutine filter_connected_i_H_psi0_SC2(key1,key2,Nint,sze,idx,idx_repeat)
endif
else
print*,'more than a double excitation, can not apply the '
print*,'SC2 dressing of the diagonal element .....'
print*,'stop !!'
print*,'degree = ',degree
stop
! print*,'more than a double excitation, can not apply the '
! print*,'SC2 dressing of the diagonal element .....'
! print*,'stop !!'
! print*,'degree = ',degree
! stop
idx(0) = 0
idx_repeat(0) = 0
endif
idx(0) = l-1
idx_repeat(0) = l_repeat-1

18
src/MonoInts/README.rst
View File

@ -91,34 +91,34 @@ Documentation
`ao_nucl_elec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L1>`_
interaction nuclear electron
`give_polynom_mult_center_mono_elec <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L157>`_
`give_polynom_mult_center_mono_elec <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L161>`_
Undocumented
`i_x1_pol_mult_mono_elec <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L285>`_
`i_x1_pol_mult_mono_elec <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L289>`_
Undocumented
`i_x2_pol_mult_mono_elec <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L357>`_
`i_x2_pol_mult_mono_elec <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L361>`_
Undocumented
`int_gaus_pol <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L428>`_
`int_gaus_pol <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L432>`_
Undocumented
`nai_pol_mult <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L82>`_
Undocumented
`v_e_n <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L409>`_
`v_e_n <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L413>`_
Undocumented
`v_phi <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L473>`_
`v_phi <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L477>`_
Undocumented
`v_r <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L457>`_
`v_r <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L461>`_
Undocumented
`v_theta <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L486>`_
`v_theta <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L490>`_
Undocumented
`wallis <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L502>`_
`wallis <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L506>`_
Undocumented
`mo_nucl_elec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_mo_ints.irp.f#L1>`_

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

@ -116,6 +116,10 @@ include 'constants.F'
enddo
const_factor = dist*rho
const = p * dist_integral
if(const_factor.ge.80.d0)then
NAI_pol_mult = 0.d0
return
endif
factor = dexp(-const_factor)
coeff = dtwo_pi * factor * p_inv
lmax = 20

2
src/NEEDED_MODULES
View File

@ -1 +1 @@
AOs Bitmask Electrons Ezfio_files MOs Nuclei Output Utils Hartree_Fock BiInts MonoInts MOGuess Dets DensityMatrix CISD Perturbation SingleRefMethod CISD_selected Selectors_full MP2 Generators_full Full_CI
AOs Bitmask Electrons Ezfio_files MOs Nuclei Output Utils Hartree_Fock BiInts MonoInts MOGuess Dets DensityMatrix CISD Perturbation SingleRefMethod CISD_selected Selectors_full MP2 Generators_full Full_CI CISD_SC2_selected CISD_SC2_selected CISD_SC2

187
src/Perturbation/epstein_nesbet.irp.f
View File

@ -24,13 +24,14 @@ subroutine pt2_epstein_nesbet(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_s
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)
do i =1,N_st
H_pert_diag(i) = h
if (dabs(CI_electronic_energy(i) - h) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - h)
H_pert_diag(i) = h*c_pert(i)*c_pert(i)
e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
else
c_pert(i) = 0.d0
e_2_pert(i) = 0.d0
c_pert(i) = -1.d0
e_2_pert(i) = -dabs(i_H_psi_array(i))
H_pert_diag(i) = h
endif
enddo
@ -62,190 +63,14 @@ subroutine pt2_epstein_nesbet_2x2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet
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)
do i =1,N_st
H_pert_diag(i) = h
delta_e = h - CI_electronic_energy(i)
e_2_pert(i) = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i)))
if (dabs(i_H_psi_array(i)) > 1.d-6) then
c_pert(i) = e_2_pert(i)/i_H_psi_array(i)
else
c_pert(i) = 0.d0
c_pert(i) = -1.d0
endif
H_pert_diag(i) = h*c_pert(i)*c_pert(i)
enddo
end
subroutine pt2_epstein_nesbet_SC2(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)
integer :: idx_repeat(0:ndet)
BEGIN_DOC
! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
!
! for the various N_st states,
!
! but with the correction in the denominator
!
! comming from the interaction of that determinant with all the others determinants
!
! that can be repeated by repeating all the double excitations
!
! : you repeat all the correlation energy already taken into account in CI_electronic_energy(1)
!
! that could be repeated to this determinant.
!
! <det_pert|H|det_pert> ---> <det_pert|H|det_pert> + delta_e_corr
!
! 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,accu_e_corr,hij,h
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array,idx_repeat)
accu_e_corr = 0.d0
do i = 1, idx_repeat(0)
call i_H_j(psi_selectors(1,1,idx_repeat(i)),det_pert,Nint,hij)
accu_e_corr = accu_e_corr + hij * psi_selectors_coef(idx_repeat(i),1)
enddo
accu_e_corr = accu_e_corr / psi_selectors_coef(1,1)
h = diag_H_mat_elem(det_pert,Nint) + accu_e_corr
do i =1,N_st
H_pert_diag(i) = h
e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
if (dabs(CI_electronic_energy(i) - h) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - h)
else
c_pert(i) = 0.d0
endif
enddo
end
subroutine pt2_epstein_nesbet_2x2_SC2(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)
integer :: idx_repeat(0:ndet)
BEGIN_DOC
! compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution
!
! for the various N_st states.
!
! but with the correction in the denominator
!
! comming from the interaction of that determinant with all the others determinants
!
! that can be repeated by repeating all the double excitations
!
! : you repeat all the correlation energy already taken into account in CI_electronic_energy(1)
!
! that could be repeated to this determinant.
!
! <det_pert|H|det_pert> ---> <det_pert|H|det_pert> + delta_e_corr
!
! e_2_pert(i) = 0.5 * (( <det_pert|H|det_pert> - E(i) ) - sqrt( ( <det_pert|H|det_pert> - E(i)) ^2 + 4 <psi(i)|H|det_pert>^2 )
!
! c_pert(i) = e_2_pert(i)/ <psi(i)|H|det_pert>
!
END_DOC
integer :: i,j
double precision :: diag_H_mat_elem,accu_e_corr,hij,delta_e,h
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array,idx_repeat)
accu_e_corr = 0.d0
do i = 1, idx_repeat(0)
call i_H_j(psi_selectors(1,1,idx_repeat(i)),det_pert,Nint,hij)
accu_e_corr = accu_e_corr + hij * psi_selectors_coef(idx_repeat(i),1)
enddo
accu_e_corr = accu_e_corr / psi_selectors_coef(1,1)
h = diag_H_mat_elem(det_pert,Nint) + accu_e_corr
do i =1,N_st
H_pert_diag(i) = h
delta_e = h - CI_electronic_energy(i)
e_2_pert(i) = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i)))
if (dabs(i_H_psi_array(i)) > 1.d-6) then
c_pert(i) = e_2_pert(i)/i_H_psi_array(i)
else
c_pert(i) = 0.d0
endif
enddo
end
subroutine pt2_epstein_nesbet_SC2_projected(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)
integer :: idx_repeat(0:ndet)
BEGIN_DOC
! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
!
! for the various N_st states,
!
! but with the correction in the denominator
!
! comming from the interaction of that determinant with all the others determinants
!
! that can be repeated by repeating all the double excitations
!
! : you repeat all the correlation energy already taken into account in CI_electronic_energy(1)
!
! that could be repeated to this determinant.
!
! BUT on the contrary with ""pt2_epstein_nesbet_SC2"", you compute the energy by projection
!
! <det_pert|H|det_pert> ---> <det_pert|H|det_pert> + delta_e_corr
!
! c_pert(1) = 1/c_HF <psi(i)|H|det_pert>/( E(i) - (<det_pert|H|det_pert> ) )
!
! e_2_pert(1) = <HF|H|det_pert> c_pert(1)
!
! NOTE :::: if you satisfy Brillouin Theorem, the singles don't contribute !!
!
END_DOC
integer :: i,j
double precision :: diag_H_mat_elem,accu_e_corr,hij,h0j,h
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array,idx_repeat)
accu_e_corr = 0.d0
call i_H_j(ref_bitmask,det_pert,Nint,h0j)
do i = 1, idx_repeat(0)
call i_H_j(psi_selectors(1,1,idx_repeat(i)),det_pert,Nint,hij)
accu_e_corr = accu_e_corr + hij * psi_selectors_coef(idx_repeat(i),1)
enddo
accu_e_corr = accu_e_corr / psi_selectors_coef(1,1)
h = diag_H_mat_elem(det_pert,Nint) + accu_e_corr
c_pert(1) = 1.d0/psi_selectors_coef(1,1) * i_H_psi_array(1) / (CI_electronic_energy(i) - h)
e_2_pert(1) = c_pert(i) * h0j
do i =2,N_st
H_pert_diag(i) = h
if (dabs(CI_electronic_energy(i) - h) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - h)
else
c_pert(i) = 0.d0
endif
e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
enddo
end

94
src/Perturbation/pert_sc2.irp.f Normal file
View File

@ -0,0 +1,94 @@
subroutine pt2_epstein_nesbet_SC2_projected(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)
integer :: idx_repeat(0:ndet)
BEGIN_DOC
! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
!
! for the various N_st states,
!
! but with the correction in the denominator
!
! comming from the interaction of that determinant with all the others determinants
!
! that can be repeated by repeating all the double excitations
!
! : you repeat all the correlation energy already taken into account in CI_electronic_energy(1)
!
! that could be repeated to this determinant.
!
! In addition, for the perturbative energetic contribution you have the standard second order
!
! e_2_pert = <psi_i|H|det_pert>^2/(Delta_E)
!
! and also the purely projected contribution
!
! H_pert_diag = <HF|H|det_pert> c_pert
END_DOC
integer :: i,j,degree
double precision :: diag_H_mat_elem,accu_e_corr,hij,h0j,h,delta_E
double precision :: repeat_all_e_corr,accu_e_corr_tmp
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array,idx_repeat)
accu_e_corr = 0.d0
call i_H_j(ref_bitmask,det_pert,Nint,h0j)
do i = 1, idx_repeat(0)
accu_e_corr = accu_e_corr + E_corr_per_selectors(idx_repeat(i))
enddo
h = diag_H_mat_elem(det_pert,Nint) + accu_e_corr
delta_E = (CI_SC2_electronic_energy(1) - h)
delta_E = 1.d0/delta_E
c_pert(1) = i_H_psi_array(1) * delta_E
e_2_pert(1) = i_H_psi_array(1) * c_pert(1)
H_pert_diag(1) = c_pert(1) * h0j/coef_hf_selector
do i =2,N_st
H_pert_diag(i) = h
if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (CI_SC2_electronic_energy(i) - h)
e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
else
c_pert(i) = -1.d0
e_2_pert(i) = -dabs(i_H_psi_array(i))
endif
enddo
end
double precision function repeat_all_e_corr(key_in)
implicit none
integer(bit_kind), intent(in) :: key_in(N_int,2)
integer :: i,degree
double precision :: accu
use bitmasks
accu = 0.d0
call get_excitation_degree(key_in,ref_bitmask,degree,N_int)
if(degree==2)then
do i = 1, N_det_selectors
call get_excitation_degree(ref_bitmask,psi_selectors(1,1,i),degree,N_int)
if(degree.ne.2)cycle
call get_excitation_degree(key_in,psi_selectors(1,1,i),degree,N_int)
if (degree<=3)cycle
accu += E_corr_per_selectors(i)
enddo
elseif(degree==1)then
do i = 1, N_det_selectors
call get_excitation_degree(ref_bitmask,psi_selectors(1,1,i),degree,N_int)
if(degree.ne.2)cycle
call get_excitation_degree(key_in,psi_selectors(1,1,i),degree,N_int)
if (degree<=2)cycle
accu += E_corr_per_selectors(i)
enddo
endif
repeat_all_e_corr = accu
end

2
src/Perturbation/perturbation_template.f
View File

@ -38,7 +38,7 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c
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) += c_pert(k) * c_pert(k) * H_pert_diag(k)
sum_H_pert_diag(k) += H_pert_diag(k)
enddo
enddo

8
src/README.rst
View File

@ -131,6 +131,14 @@ Needed Modules
* `CISD <http://github.com/LCPQ/quantum_package/tree/master/src/CISD>`_
* `Perturbation <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation>`_
* `SingleRefMethod <http://github.com/LCPQ/quantum_package/tree/master/src/SingleRefMethod>`_
* `CISD_selected <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_selected>`_
* `Selectors_full <http://github.com/LCPQ/quantum_package/tree/master/src/Selectors_full>`_
* `MP2 <http://github.com/LCPQ/quantum_package/tree/master/src/MP2>`_
* `Generators_full <http://github.com/LCPQ/quantum_package/tree/master/src/Generators_full>`_
* `Full_CI <http://github.com/LCPQ/quantum_package/tree/master/src/Full_CI>`_
* `CISD_SC2_selected <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_SC2_selected>`_
* `CISD_SC2_selected <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_SC2_selected>`_
* `CISD_SC2 <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_SC2>`_
Documentation
=============

59
src/Selectors_full/e_corr_selectors.irp.f Normal file
View File

@ -0,0 +1,59 @@
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, E_corr_per_selectors, (N_det_selectors)]
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
integer :: i,degree
double precision :: hij,inv_selectors_coef_hf
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)
E_corr_per_selectors(i) = psi_selectors_coef(i,1) * hij
elseif(exc_degree_per_selectors(i) == 0)then
coef_hf_selector = psi_selectors_coef(i,1)
E_corr_per_selectors(i) = -1000.d0
else
E_corr_per_selectors(i) = -1000.d0
endif
enddo
inv_selectors_coef_hf = 1.d0/coef_hf_selector
do i = 1,n_double_selectors
E_corr_per_selectors(double_index_selectors(i)) *=inv_selectors_coef_hf
enddo
END_PROVIDER

20
src/Utils/README.rst
View File

@ -46,6 +46,18 @@ Documentation
m : Coefficients matrix is MxN, ( array is (LDC,N) )
.br
`abort_all <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/abort.irp.f#L1>`_
If True, all the calculation is aborted
`abort_here <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/abort.irp.f#L10>`_
If True, all the calculation is aborted
`control_c <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/abort.irp.f#L32>`_
What to do on Ctrl-C. If two Ctrl-C are pressed within 1 sec, the calculation if aborted.
`trap_signals <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/abort.irp.f#L18>`_
What to do when a signal is caught. Here, trap Ctrl-C and call the control_C subroutine.
`add_poly <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/integration.irp.f#L243>`_
Add two polynomials
D(t) =! D(t) +( B(t)+C(t))
@ -80,7 +92,7 @@ Documentation
into
fact_k (x-x_P)^iorder(1) (y-y_P)^iorder(2) (z-z_P)^iorder(3) exp(-p(r-P)^2)
`hermite <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/integration.irp.f#L468>`_
`hermite <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/integration.irp.f#L472>`_
Hermite polynomial
`multiply_poly <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/integration.irp.f#L201>`_
@ -96,13 +108,13 @@ Documentation
\int_0^1 dx \exp(-p x^2) x^n
.br
`rint1 <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/integration.irp.f#L524>`_
`rint1 <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/integration.irp.f#L528>`_
Standard version of rint
`rint_large_n <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/integration.irp.f#L493>`_
`rint_large_n <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/integration.irp.f#L497>`_
Version of rint for large values of n
`rint_sum <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/integration.irp.f#L417>`_
`rint_sum <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/integration.irp.f#L421>`_
Needed for the calculation of two-electron integrals.
`overlap_gaussian_x <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/one_e_integration.irp.f#L1>`_

15
src/Utils/integration.irp.f
View File

@ -151,7 +151,7 @@ subroutine gaussian_product(a,xa,b,xb,k,p,xp)
k=0.d0
return
endif
k = exp(-k)
k = dexp(-k)
xp(1) = (a*xa(1)+b*xb(1))*p_inv
xp(2) = (a*xa(2)+b*xb(2))*p_inv
xp(3) = (a*xa(3)+b*xb(3))*p_inv
@ -398,7 +398,11 @@ double precision function rint(n,rho)
else
if(n.le.20)then
u_inv=1.d0/dsqrt(rho)
v=dexp(-rho)
if(rho.gt.80.d0)then
v=0.d0
else
v=dexp(-rho)
endif
u=rho*u_inv
two_rho_inv = 0.5d0*u_inv*u_inv
val0=0.5d0*u_inv*sqpi*erf(u)
@ -444,7 +448,12 @@ double precision function rint_sum(n_pt_out,rho,d1)
else
v=dexp(-rho)
if(rho.gt.80.d0)then
v=0.d0
else
v=dexp(-rho)
endif
u_inv=1.d0/dsqrt(rho)
u=rho*u_inv
two_rho_inv = 0.5d0*u_inv*u_inv

6
src/Utils/one_e_integration.irp.f
View File

@ -55,7 +55,7 @@ subroutine overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,&
integer :: iorder_p(3)
call give_explicit_poly_and_gaussian(P_new,P_center,p,fact_p,iorder_p,alpha,beta,power_A,power_B,A_center,B_center,dim)
if(fact_p.lt.0.000001d0)then
if(fact_p.lt.1d-10)then
overlap_x = 0.d0
overlap_y = 0.d0
overlap_z = 0.d0
@ -122,6 +122,10 @@ subroutine overlap_x_abs(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,
rho = alpha * beta * p_inv
dist = (A_center - B_center)*(A_center - B_center)
P_center = (alpha * A_center + beta * B_center) * p_inv
if(rho*dist.gt.80.d0)then
overlap_x= 0.d0
return
endif
factor = dexp(-rho * dist)
double precision :: tmp