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

Merge branch 'master' of github.com:LCPQ/quantum_package

Browse Source 
Conflicts:
	src/Dets/README.rst
	src/NEEDED_MODULES
This commit is contained in:
Anthony Scemama 2014-05-30 23:20:52 +02:00
commit 3360bdf993
29 changed files with 472 additions and 295 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
scripts/generate_h_apply.py
View File

@ -124,6 +124,7 @@ class H_apply(object):
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

7
src/TestFock/README.rst → src/CISD_SC2/README.rst
View File

@ -1,8 +1,7 @@
===============
TestFock Module
CISD_SC2 Module
===============
* This is a test
Documentation
=============
@ -20,6 +19,8 @@ Needed Modules
* `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>`_
@ -27,5 +28,7 @@ Needed Modules
* `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

28
src/Dets/README.rst
View File

@ -77,12 +77,6 @@ Documentation
`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
`det_connections <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/connections.irp.f#L10>`_
.br
`n_con_int <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/connections.irp.f#L2>`_
Number of integers to represent the connections between determinants
`davidson_converged <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L383>`_
True if the Davidson algorithm is converged
@ -213,7 +207,7 @@ Documentation
.br
idx(0) is the number of determinants that interact with key1
`filter_connected_davidson <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/filter_connected.irp.f#L95>`_
`filter_connected_davidson <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/filter_connected.irp.f#L101>`_
Filters out the determinants that are not connected by H
.br
returns the array idx which contains the index of the
@ -224,7 +218,7 @@ Documentation
.br
idx(0) is the number of determinants that interact with key1
`filter_connected_i_h_psi0 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/filter_connected.irp.f#L211>`_
`filter_connected_i_h_psi0 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/filter_connected.irp.f#L233>`_
returns the array idx which contains the index of the
.br
determinants in the array key1 that interact
@ -233,8 +227,16 @@ Documentation
.br
idx(0) is the number of determinants that interact with key1
`filter_connected_i_h_psi0_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/filter_connected.irp.f#L310>`_
Undocumented
`filter_connected_i_h_psi0_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/filter_connected.irp.f#L332>`_
standard filter_connected_i_H_psi but returns in addition
.br
the array of the index of the non connected determinants to key1
.br
in order to know what double excitation can be repeated on key1
.br
idx_repeat(0) is the number of determinants that can be used
.br
to repeat the excitations
`get_s2 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/s2.irp.f#L1>`_
Returns <S^2>
@ -261,6 +263,9 @@ Documentation
s1,s2 : Spins (1:alpha, 2:beta)
degree : Degree of excitation
`det_connections <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L898>`_
.br
`diag_h_mat_elem <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L659>`_
Computes <i|H|i>
@ -308,6 +313,9 @@ Documentation
.br
to repeat the excitations
`n_con_int <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L890>`_
Number of integers to represent the connections between determinants
`h_matrix_all_dets <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/utils.irp.f#L1>`_
H matrix on the basis of the slater deter;inants defined by psi_det

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

@ -360,6 +360,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
@ -374,7 +375,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
@ -393,7 +394,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
@ -414,7 +415,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
@ -438,7 +439,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
@ -531,11 +532,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 TestFock
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
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)
H_pert_diag(i) = h*c_pert(i)*c_pert(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>`_

11
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)
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
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