From 868ef807bd654d7742a01434a51ff36a37231288 Mon Sep 17 00:00:00 2001 From: Manu Date: Fri, 30 May 2014 18:07:04 +0200 Subject: [PATCH] CISD_SC2 ok, CISD_SC2_selected ok --- scripts/generate_h_apply.py | 5 +- src/CISD/README.rst | 21 +- src/CISD/cisd_lapack.irp.f | 7 +- src/CISD/cisd_sc2.irp.f | 23 --- src/CISD_SC2/README.rst | 34 ++++ src/{CISD => CISD_SC2}/SC2.irp.f | 63 ++++-- src/CISD_SC2/cisd_SC2.irp.f | 14 ++ src/CISD_SC2/diagonalize_CI_SC2.irp.f | 69 +++++++ src/CISD_SC2_selected/ASSUMPTIONS.rst | 0 src/CISD_SC2_selected/Makefile | 8 + src/CISD_SC2_selected/NEEDED_MODULES | 1 + src/CISD_SC2_selected/README.rst | 37 ++++ .../cisd_sc2_selection.irp.f | 32 +++ src/CISD_selected/NEEDED_MODULES | 2 +- src/CISD_selected/README.rst | 1 + src/CISD_selected/cisd_selection.irp.f | 13 +- src/Dets/Makefile | 4 - src/Dets/README.rst | 64 ++++-- src/Dets/filter_connected.irp.f | 21 +- src/MonoInts/README.rst | 18 +- src/MonoInts/pot_ao_ints.irp.f | 4 + src/NEEDED_MODULES | 2 +- src/Perturbation/epstein_nesbet.irp.f | 187 +----------------- src/Perturbation/pert_sc2.irp.f | 94 +++++++++ src/Perturbation/perturbation_template.f | 2 +- src/README.rst | 8 + src/Selectors_full/e_corr_selectors.irp.f | 59 ++++++ src/Utils/README.rst | 20 +- src/Utils/integration.irp.f | 15 +- src/Utils/one_e_integration.irp.f | 6 +- 30 files changed, 532 insertions(+), 302 deletions(-) delete mode 100644 src/CISD/cisd_sc2.irp.f create mode 100644 src/CISD_SC2/README.rst rename src/{CISD => CISD_SC2}/SC2.irp.f (74%) create mode 100644 src/CISD_SC2/cisd_SC2.irp.f create mode 100644 src/CISD_SC2/diagonalize_CI_SC2.irp.f create mode 100644 src/CISD_SC2_selected/ASSUMPTIONS.rst create mode 100644 src/CISD_SC2_selected/Makefile create mode 100644 src/CISD_SC2_selected/NEEDED_MODULES create mode 100644 src/CISD_SC2_selected/README.rst create mode 100644 src/CISD_SC2_selected/cisd_sc2_selection.irp.f create mode 100644 src/Perturbation/pert_sc2.irp.f create mode 100644 src/Selectors_full/e_corr_selectors.irp.f diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index eeac4669..a6143b23 100755 --- a/scripts/generate_h_apply.py +++ b/scripts/generate_h_apply.py @@ -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, & diff --git a/src/CISD/README.rst b/src/CISD/README.rst index 5d4c9b4d..71c0f188 100644 --- a/src/CISD/README.rst +++ b/src/CISD/README.rst @@ -35,26 +35,7 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`cisd_sc2 `_ - 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 `_ - Undocumented - -`cisd `_ +`cisd `_ Undocumented diff --git a/src/CISD/cisd_lapack.irp.f b/src/CISD/cisd_lapack.irp.f index c57d11f1..374dc9c2 100644 --- a/src/CISD/cisd_lapack.irp.f +++ b/src/CISD/cisd_lapack.irp.f @@ -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 diff --git a/src/CISD/cisd_sc2.irp.f b/src/CISD/cisd_sc2.irp.f deleted file mode 100644 index 6e4b1bd3..00000000 --- a/src/CISD/cisd_sc2.irp.f +++ /dev/null @@ -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 diff --git a/src/CISD_SC2/README.rst b/src/CISD_SC2/README.rst new file mode 100644 index 00000000..938b01aa --- /dev/null +++ b/src/CISD_SC2/README.rst @@ -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 `_ +* `BiInts `_ +* `Bitmask `_ +* `CISD `_ +* `Dets `_ +* `Electrons `_ +* `Ezfio_files `_ +* `Hartree_Fock `_ +* `MonoInts `_ +* `MOs `_ +* `Nuclei `_ +* `Output `_ +* `SingleRefMethod `_ +* `Utils `_ +* `Selectors_full `_ + diff --git a/src/CISD/SC2.irp.f b/src/CISD_SC2/SC2.irp.f similarity index 74% rename from src/CISD/SC2.irp.f rename to src/CISD_SC2/SC2.irp.f index 35f19e2b..64eb3934 100644 --- a/src/CISD/SC2.irp.f +++ b/src/CISD_SC2/SC2.irp.f @@ -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 + diff --git a/src/CISD_SC2/cisd_SC2.irp.f b/src/CISD_SC2/cisd_SC2.irp.f new file mode 100644 index 00000000..6cae2e9c --- /dev/null +++ b/src/CISD_SC2/cisd_SC2.irp.f @@ -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 diff --git a/src/CISD_SC2/diagonalize_CI_SC2.irp.f b/src/CISD_SC2/diagonalize_CI_SC2.irp.f new file mode 100644 index 00000000..cd01296a --- /dev/null +++ b/src/CISD_SC2/diagonalize_CI_SC2.irp.f @@ -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 diff --git a/src/CISD_SC2_selected/ASSUMPTIONS.rst b/src/CISD_SC2_selected/ASSUMPTIONS.rst new file mode 100644 index 00000000..e69de29b diff --git a/src/CISD_SC2_selected/Makefile b/src/CISD_SC2_selected/Makefile new file mode 100644 index 00000000..b2ea1de1 --- /dev/null +++ b/src/CISD_SC2_selected/Makefile @@ -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 diff --git a/src/CISD_SC2_selected/NEEDED_MODULES b/src/CISD_SC2_selected/NEEDED_MODULES new file mode 100644 index 00000000..0b32345c --- /dev/null +++ b/src/CISD_SC2_selected/NEEDED_MODULES @@ -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 diff --git a/src/CISD_SC2_selected/README.rst b/src/CISD_SC2_selected/README.rst new file mode 100644 index 00000000..e815aa1b --- /dev/null +++ b/src/CISD_SC2_selected/README.rst @@ -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 `_ +* `BiInts `_ +* `Bitmask `_ +* `CISD `_ +* `CISD_SC2 `_ +* `CISD_selected `_ +* `Dets `_ +* `Electrons `_ +* `Ezfio_files `_ +* `Hartree_Fock `_ +* `MonoInts `_ +* `MOs `_ +* `Nuclei `_ +* `Output `_ +* `Perturbation `_ +* `Selectors_full `_ +* `SingleRefMethod `_ +* `Utils `_ + diff --git a/src/CISD_SC2_selected/cisd_sc2_selection.irp.f b/src/CISD_SC2_selected/cisd_sc2_selection.irp.f new file mode 100644 index 00000000..96d17199 --- /dev/null +++ b/src/CISD_SC2_selected/cisd_sc2_selection.irp.f @@ -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 diff --git a/src/CISD_selected/NEEDED_MODULES b/src/CISD_selected/NEEDED_MODULES index ac42ef59..48b7f863 100644 --- a/src/CISD_selected/NEEDED_MODULES +++ b/src/CISD_selected/NEEDED_MODULES @@ -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 diff --git a/src/CISD_selected/README.rst b/src/CISD_selected/README.rst index 4fed26e5..eee8e392 100644 --- a/src/CISD_selected/README.rst +++ b/src/CISD_selected/README.rst @@ -26,6 +26,7 @@ Needed Modules * `BiInts `_ * `Bitmask `_ * `CISD `_ +* `CISD_SC2 `_ * `Dets `_ * `Electrons `_ * `Ezfio_files `_ diff --git a/src/CISD_selected/cisd_selection.irp.f b/src/CISD_selected/cisd_selection.irp.f index 6d88f61e..da73692e 100644 --- a/src/CISD_selected/cisd_selection.irp.f +++ b/src/CISD_selected/cisd_selection.irp.f @@ -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 diff --git a/src/Dets/Makefile b/src/Dets/Makefile index ea6bc62c..f97f77c0 100644 --- a/src/Dets/Makefile +++ b/src/Dets/Makefile @@ -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 diff --git a/src/Dets/README.rst b/src/Dets/README.rst index 7c7ffade..54dfd5b2 100644 --- a/src/Dets/README.rst +++ b/src/Dets/README.rst @@ -51,9 +51,10 @@ Documentation .. NEEDED_MODULES file. `copy_h_apply_buffer_to_wf `_ - 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 `_ +`fill_h_apply_buffer_no_selection `_ Fill the H_apply buffer with determiants for CISD `h_apply_buffer_allocated `_ @@ -69,18 +70,18 @@ Documentation `connected_to_ref `_ Undocumented -`det_is_not_or_may_be_in_ref `_ +`det_is_not_or_may_be_in_ref `_ If true, det is not in ref If false, det may be in ref -`key_pattern_not_in_ref `_ +`key_pattern_not_in_ref `_ Min and max values of the integers of the keys of the reference -`davidson_converged `_ +`davidson_converged `_ True if the Davidson algorithm is converged -`davidson_criterion `_ - Can be : [ energy | residual | both ] +`davidson_criterion `_ + Can be : [ energy | residual | both | wall_time | cpu_time | iterations ] `davidson_diag `_ Davidson diagonalization. @@ -126,24 +127,41 @@ Documentation `davidson_sze_max `_ Max number of Davidson sizes -`davidson_threshold `_ - Can be : [ energy | residual | both ] +`davidson_threshold `_ + Can be : [ energy | residual | both | wall_time | cpu_time | iterations ] -`n_det `_ +`n_det `_ Number of determinants in the wave function +`n_det_reference `_ + Number of determinants in the reference wave function + `n_states `_ Number of states to consider -`psi_coef `_ - The wave function. Initialized with Hartree-Fock +`psi_average_norm_contrib `_ + Contribution of determinants to the state-averaged density -`psi_det `_ - The wave function. Initialized with Hartree-Fock +`psi_average_norm_contrib_sorted `_ + Wave function sorted by determinants (state-averaged) -`psi_det_size `_ +`psi_coef `_ + The wave function. Initialized with Hartree-Fock if the EZFIO file + is empty + +`psi_coef_sorted `_ + Wave function sorted by determinants (state-averaged) + +`psi_det `_ + The wave function. Initialized with Hartree-Fock if the EZFIO file + is empty + +`psi_det_size `_ Size of the psi_det/psi_coef arrays +`psi_det_sorted `_ + Wave function sorted by determinants (state-averaged) + `double_exc_bitmask `_ 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 `_ + Eigenvectors/values of the CI matrix + +`ci_electronic_energy `_ + Eigenvectors/values of the CI matrix + +`ci_energy `_ + N_states lowest eigenvalues of the CI matrix + +`diag_algorithm `_ + Diagonalization algorithm (Davidson or Lapack) + +`diagonalize_ci `_ + Replace the coefficients of the CI states by the coefficients of the + eigenstates of the CI matrix + `filter_connected `_ Filters out the determinants that are not connected by H .br diff --git a/src/Dets/filter_connected.irp.f b/src/Dets/filter_connected.irp.f index c9814009..3bda30e6 100644 --- a/src/Dets/filter_connected.irp.f +++ b/src/Dets/filter_connected.irp.f @@ -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 diff --git a/src/MonoInts/README.rst b/src/MonoInts/README.rst index 052e2053..a9bd07bc 100644 --- a/src/MonoInts/README.rst +++ b/src/MonoInts/README.rst @@ -91,34 +91,34 @@ Documentation `ao_nucl_elec_integral `_ interaction nuclear electron -`give_polynom_mult_center_mono_elec `_ +`give_polynom_mult_center_mono_elec `_ Undocumented -`i_x1_pol_mult_mono_elec `_ +`i_x1_pol_mult_mono_elec `_ Undocumented -`i_x2_pol_mult_mono_elec `_ +`i_x2_pol_mult_mono_elec `_ Undocumented -`int_gaus_pol `_ +`int_gaus_pol `_ Undocumented `nai_pol_mult `_ Undocumented -`v_e_n `_ +`v_e_n `_ Undocumented -`v_phi `_ +`v_phi `_ Undocumented -`v_r `_ +`v_r `_ Undocumented -`v_theta `_ +`v_theta `_ Undocumented -`wallis `_ +`wallis `_ Undocumented `mo_nucl_elec_integral `_ diff --git a/src/MonoInts/pot_ao_ints.irp.f b/src/MonoInts/pot_ao_ints.irp.f index 3aedeebc..deac43c7 100644 --- a/src/MonoInts/pot_ao_ints.irp.f +++ b/src/MonoInts/pot_ao_ints.irp.f @@ -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 diff --git a/src/NEEDED_MODULES b/src/NEEDED_MODULES index 6151539c..c1fcfc03 100644 --- a/src/NEEDED_MODULES +++ b/src/NEEDED_MODULES @@ -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 diff --git a/src/Perturbation/epstein_nesbet.irp.f b/src/Perturbation/epstein_nesbet.irp.f index d6a903c7..36d9ab1f 100644 --- a/src/Perturbation/epstein_nesbet.irp.f +++ b/src/Perturbation/epstein_nesbet.irp.f @@ -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. - ! - ! ---> + delta_e_corr - ! - ! c_pert(i) = /( E(i) - ( ) ) - ! - ! e_2_pert(i) = ^2/( E(i) - ( ) ) - ! - 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. - ! - ! ---> + delta_e_corr - ! - ! e_2_pert(i) = 0.5 * (( - E(i) ) - sqrt( ( - E(i)) ^2 + 4 ^2 ) - ! - ! c_pert(i) = e_2_pert(i)/ - ! - 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 - ! - ! ---> + delta_e_corr - ! - ! c_pert(1) = 1/c_HF /( E(i) - ( ) ) - ! - ! e_2_pert(1) = 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 diff --git a/src/Perturbation/pert_sc2.irp.f b/src/Perturbation/pert_sc2.irp.f new file mode 100644 index 00000000..d9c99cd2 --- /dev/null +++ b/src/Perturbation/pert_sc2.irp.f @@ -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 = ^2/(Delta_E) + ! + ! and also the purely projected contribution + ! + ! H_pert_diag = 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 diff --git a/src/Perturbation/perturbation_template.f b/src/Perturbation/perturbation_template.f index e3033820..d96b557a 100644 --- a/src/Perturbation/perturbation_template.f +++ b/src/Perturbation/perturbation_template.f @@ -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 diff --git a/src/README.rst b/src/README.rst index 5882eb49..4b6fab44 100644 --- a/src/README.rst +++ b/src/README.rst @@ -131,6 +131,14 @@ Needed Modules * `CISD `_ * `Perturbation `_ * `SingleRefMethod `_ +* `CISD_selected `_ +* `Selectors_full `_ +* `MP2 `_ +* `Generators_full `_ +* `Full_CI `_ +* `CISD_SC2_selected `_ +* `CISD_SC2_selected `_ +* `CISD_SC2 `_ Documentation ============= diff --git a/src/Selectors_full/e_corr_selectors.irp.f b/src/Selectors_full/e_corr_selectors.irp.f new file mode 100644 index 00000000..580c478b --- /dev/null +++ b/src/Selectors_full/e_corr_selectors.irp.f @@ -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) = * 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 diff --git a/src/Utils/README.rst b/src/Utils/README.rst index a93db8dd..3fd5bd8d 100644 --- a/src/Utils/README.rst +++ b/src/Utils/README.rst @@ -46,6 +46,18 @@ Documentation m : Coefficients matrix is MxN, ( array is (LDC,N) ) .br +`abort_all `_ + If True, all the calculation is aborted + +`abort_here `_ + If True, all the calculation is aborted + +`control_c `_ + What to do on Ctrl-C. If two Ctrl-C are pressed within 1 sec, the calculation if aborted. + +`trap_signals `_ + What to do when a signal is caught. Here, trap Ctrl-C and call the control_C subroutine. + `add_poly `_ 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 `_ +`hermite `_ Hermite polynomial `multiply_poly `_ @@ -96,13 +108,13 @@ Documentation \int_0^1 dx \exp(-p x^2) x^n .br -`rint1 `_ +`rint1 `_ Standard version of rint -`rint_large_n `_ +`rint_large_n `_ Version of rint for large values of n -`rint_sum `_ +`rint_sum `_ Needed for the calculation of two-electron integrals. `overlap_gaussian_x `_ diff --git a/src/Utils/integration.irp.f b/src/Utils/integration.irp.f index 2f06687d..e1a6a59e 100644 --- a/src/Utils/integration.irp.f +++ b/src/Utils/integration.irp.f @@ -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 diff --git a/src/Utils/one_e_integration.irp.f b/src/Utils/one_e_integration.irp.f index f3a7e0a6..d7e52d8e 100644 --- a/src/Utils/one_e_integration.irp.f +++ b/src/Utils/one_e_integration.irp.f @@ -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