From 503cea9378c483b79214e67ffabe493fb7531d0a Mon Sep 17 00:00:00 2001 From: Manu Date: Fri, 6 Jun 2014 16:22:54 +0200 Subject: [PATCH 1/4] add threshold_convergence_SC2 --- src/CISD_SC2_selected/cisd_sc2_selection.irp.f | 6 +++--- src/Dets/diagonalize_CI_SC2.irp.f | 10 +++++++++- 2 files changed, 12 insertions(+), 4 deletions(-) diff --git a/src/CISD_SC2_selected/cisd_sc2_selection.irp.f b/src/CISD_SC2_selected/cisd_sc2_selection.irp.f index 1bb9a389..3610a0fe 100644 --- a/src/CISD_SC2_selected/cisd_sc2_selection.irp.f +++ b/src/CISD_SC2_selected/cisd_sc2_selection.irp.f @@ -13,9 +13,9 @@ program cisd_sc2_selected pt2 = 1.d0 perturbation = "epstein_nesbet_sc2_projected" E_old(1) = HF_energy - davidson_threshold = 1.d-4 + davidson_threshold = 1.d-6 - do while (maxval(abs(pt2(1:N_st))) > 1.d-6) + do while (maxval(abs(pt2(1:N_st))) > 1.d-4) print*,'----' print*,'' call H_apply_cisd_selection(perturbation,pt2, norm_pert, H_pert_diag, N_st) @@ -39,7 +39,7 @@ program cisd_sc2_selected enddo pt2 = 0.d0 call H_apply_PT2(pt2, norm_pert, H_pert_diag, N_st) - davidson_threshold = 1.d-8 + davidson_threshold = 1.d-10 touch davidson_threshold davidson_criterion do i = 1, N_st max = 0.d0 diff --git a/src/Dets/diagonalize_CI_SC2.irp.f b/src/Dets/diagonalize_CI_SC2.irp.f index 6873f930..c1b67e54 100644 --- a/src/Dets/diagonalize_CI_SC2.irp.f +++ b/src/Dets/diagonalize_CI_SC2.irp.f @@ -15,6 +15,14 @@ BEGIN_PROVIDER [ double precision, CI_SC2_energy, (N_states) ] END_PROVIDER + BEGIN_PROVIDER [ double precision, threshold_convergence_SC2] + implicit none + BEGIN_DOC + ! convergence of the correlation energy of SC2 iterations + END_DOC + threshold_convergence_SC2 = 1.d-8 + + 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 @@ -32,7 +40,7 @@ END_PROVIDER double precision :: convergence call CISD_SC2(psi_det,CI_SC2_eigenvectors,CI_SC2_electronic_energy, & - size(CI_SC2_eigenvectors,1),N_det,N_states,N_int,davidson_threshold) + size(CI_SC2_eigenvectors,1),N_det,N_states,N_int,threshold_convergence_SC2) END_PROVIDER subroutine diagonalize_CI_SC2 From 09b43e57e0d1066a4e1e5c5ebaef2c80f9f0dad1 Mon Sep 17 00:00:00 2001 From: Manu Date: Wed, 11 Jun 2014 00:44:07 +0200 Subject: [PATCH 2/4] add CID, CID_selected, CID_SC2_selected, just for fun ... --- src/CID/ASSUMPTIONS.rst | 0 src/CID/H_apply.irp.f | 9 +++ src/CID/Makefile | 8 +++ src/CID/NEEDED_MODULES | 2 + src/CID/README.rst | 43 ++++++++++++ src/CID/cid.irp.f | 20 ++++++ src/CID/cid_lapack.irp.f | 16 +++++ src/CID/do_mono_double.irp.f | 19 +++++ src/CID_SC2_selected/ASSUMPTIONS.rst | 0 src/CID_SC2_selected/H_apply.irp.f | 10 +++ src/CID_SC2_selected/Makefile | 8 +++ src/CID_SC2_selected/NEEDED_MODULES | 1 + src/CID_SC2_selected/README.rst | 40 +++++++++++ src/CID_SC2_selected/cid_sc2_selection.irp.f | 73 ++++++++++++++++++++ src/CID_SC2_selected/do_mono_double.irp.f | 19 +++++ src/CID_selected/ASSUMPTIONS.rst | 0 src/CID_selected/H_apply.irp.f | 31 +++++++++ src/CID_selected/Makefile | 8 +++ src/CID_selected/NEEDED_MODULES | 1 + src/CID_selected/README.rst | 41 +++++++++++ src/CID_selected/cid_selection.irp.f | 34 +++++++++ src/CID_selected/do_mono_double.irp.f | 19 +++++ src/CISD/do_mono_double.irp.f | 19 +++++ src/CISD_SC2_selected/do_mono_double.irp.f | 19 +++++ src/CISD_selected/do_mono_double.irp.f | 19 +++++ src/Dets/H_apply_template.f | 22 ++++-- src/Dets/filter_connected.irp.f | 8 +-- 27 files changed, 478 insertions(+), 11 deletions(-) create mode 100644 src/CID/ASSUMPTIONS.rst create mode 100644 src/CID/H_apply.irp.f create mode 100644 src/CID/Makefile create mode 100644 src/CID/NEEDED_MODULES create mode 100644 src/CID/README.rst create mode 100644 src/CID/cid.irp.f create mode 100644 src/CID/cid_lapack.irp.f create mode 100644 src/CID/do_mono_double.irp.f create mode 100644 src/CID_SC2_selected/ASSUMPTIONS.rst create mode 100644 src/CID_SC2_selected/H_apply.irp.f create mode 100644 src/CID_SC2_selected/Makefile create mode 100644 src/CID_SC2_selected/NEEDED_MODULES create mode 100644 src/CID_SC2_selected/README.rst create mode 100644 src/CID_SC2_selected/cid_sc2_selection.irp.f create mode 100644 src/CID_SC2_selected/do_mono_double.irp.f create mode 100644 src/CID_selected/ASSUMPTIONS.rst create mode 100644 src/CID_selected/H_apply.irp.f create mode 100644 src/CID_selected/Makefile create mode 100644 src/CID_selected/NEEDED_MODULES create mode 100644 src/CID_selected/README.rst create mode 100644 src/CID_selected/cid_selection.irp.f create mode 100644 src/CID_selected/do_mono_double.irp.f create mode 100644 src/CISD/do_mono_double.irp.f create mode 100644 src/CISD_SC2_selected/do_mono_double.irp.f create mode 100644 src/CISD_selected/do_mono_double.irp.f diff --git a/src/CID/ASSUMPTIONS.rst b/src/CID/ASSUMPTIONS.rst new file mode 100644 index 00000000..e69de29b diff --git a/src/CID/H_apply.irp.f b/src/CID/H_apply.irp.f new file mode 100644 index 00000000..0df1da38 --- /dev/null +++ b/src/CID/H_apply.irp.f @@ -0,0 +1,9 @@ +! Generates subroutine H_apply_cisd +! ---------------------------------- + +BEGIN_SHELL [ /usr/bin/env python ] +from generate_h_apply import H_apply +H = H_apply("cisd") +print H +END_SHELL + diff --git a/src/CID/Makefile b/src/CID/Makefile new file mode 100644 index 00000000..b2ea1de1 --- /dev/null +++ b/src/CID/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/CID/NEEDED_MODULES b/src/CID/NEEDED_MODULES new file mode 100644 index 00000000..4b91f009 --- /dev/null +++ b/src/CID/NEEDED_MODULES @@ -0,0 +1,2 @@ +AOs BiInts Bitmask Dets Electrons Ezfio_files Hartree_Fock MonoInts MOs Nuclei Output SingleRefMethod Utils Selectors_full + diff --git a/src/CID/README.rst b/src/CID/README.rst new file mode 100644 index 00000000..003b51fa --- /dev/null +++ b/src/CID/README.rst @@ -0,0 +1,43 @@ +CISD +==== + +This is a test directory which builds a CISD by setting the follwoing rules: + +* The only generator determinant is the Hartee-Fock (single-reference method) +* All generated determinants are included in the wave function (no perturbative + selection) + +These rules are set in the ``H_apply.irp.f`` file. + +Needed Modules +============== + +.. Do not edit this section. It was auto-generated from the +.. NEEDED_MODULES file. + +* `AOs `_ +* `BiInts `_ +* `Bitmask `_ +* `Dets `_ +* `Electrons `_ +* `Ezfio_files `_ +* `Hartree_Fock `_ +* `MonoInts `_ +* `MOs `_ +* `Nuclei `_ +* `Output `_ +* `SingleRefMethod `_ +* `Utils `_ +* `Selectors_full `_ + +Documentation +============= + +.. Do not edit this section. It was auto-generated from the +.. NEEDED_MODULES file. + +`cisd `_ + Undocumented + + + diff --git a/src/CID/cid.irp.f b/src/CID/cid.irp.f new file mode 100644 index 00000000..0b6bc9fd --- /dev/null +++ b/src/CID/cid.irp.f @@ -0,0 +1,20 @@ +program cisd + implicit none + integer :: i + + print *, 'HF = ', HF_energy + print *, 'N_states = ', N_states + N_det = 1 + touch psi_det psi_coef N_det + 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 + 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/CID/cid_lapack.irp.f b/src/CID/cid_lapack.irp.f new file mode 100644 index 00000000..374dc9c2 --- /dev/null +++ b/src/CID/cid_lapack.irp.f @@ -0,0 +1,16 @@ +program cisd + implicit none + integer :: i + + diag_algorithm = "Lapack" + touch diag_algorithm + 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 + enddo + +end diff --git a/src/CID/do_mono_double.irp.f b/src/CID/do_mono_double.irp.f new file mode 100644 index 00000000..f211879e --- /dev/null +++ b/src/CID/do_mono_double.irp.f @@ -0,0 +1,19 @@ +BEGIN_PROVIDER [logical, do_double_excitations] + implicit none + BEGIN_DOC + ! if True then the double excitations are performed in the calculation + ! always true in the CISD + END_DOC + do_double_excitations = .True. + +END_PROVIDER + +BEGIN_PROVIDER [logical, do_mono_excitations] + implicit none + BEGIN_DOC + ! if True then the mono excitations are performed in the calculation + ! always true in the CISD + END_DOC + do_mono_excitations = .False. + +END_PROVIDER diff --git a/src/CID_SC2_selected/ASSUMPTIONS.rst b/src/CID_SC2_selected/ASSUMPTIONS.rst new file mode 100644 index 00000000..e69de29b diff --git a/src/CID_SC2_selected/H_apply.irp.f b/src/CID_SC2_selected/H_apply.irp.f new file mode 100644 index 00000000..79668af7 --- /dev/null +++ b/src/CID_SC2_selected/H_apply.irp.f @@ -0,0 +1,10 @@ +use bitmasks +BEGIN_SHELL [ /usr/bin/env python ] +from generate_h_apply import * +from perturbation import perturbations + +s = H_apply("PT2",SingleRef=True) +s.set_perturbation("epstein_nesbet_sc2_projected") +print s +END_SHELL + diff --git a/src/CID_SC2_selected/Makefile b/src/CID_SC2_selected/Makefile new file mode 100644 index 00000000..b2ea1de1 --- /dev/null +++ b/src/CID_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/CID_SC2_selected/NEEDED_MODULES b/src/CID_SC2_selected/NEEDED_MODULES new file mode 100644 index 00000000..bb98e702 --- /dev/null +++ b/src/CID_SC2_selected/NEEDED_MODULES @@ -0,0 +1 @@ +AOs BiInts Bitmask CISD SC2 CISD_selected Dets Electrons Ezfio_files Hartree_Fock MonoInts MOs Nuclei Output Perturbation Selectors_full SingleRefMethod Utils diff --git a/src/CID_SC2_selected/README.rst b/src/CID_SC2_selected/README.rst new file mode 100644 index 00000000..b6206850 --- /dev/null +++ b/src/CID_SC2_selected/README.rst @@ -0,0 +1,40 @@ +======================== +CISD_SC2_selected Module +======================== + +Documentation +============= + +.. Do not edit this section. It was auto-generated from the +.. NEEDED_MODULES file. + +`cisd_sc2_selected `_ + Undocumented + + + +Needed Modules +============== + +.. Do not edit this section. It was auto-generated from the +.. NEEDED_MODULES file. + +* `AOs `_ +* `BiInts `_ +* `Bitmask `_ +* `CISD `_ +* `SC2 `_ +* `CISD_selected `_ +* `Dets `_ +* `Electrons `_ +* `Ezfio_files `_ +* `Hartree_Fock `_ +* `MonoInts `_ +* `MOs `_ +* `Nuclei `_ +* `Output `_ +* `Perturbation `_ +* `Selectors_full `_ +* `SingleRefMethod `_ +* `Utils `_ + diff --git a/src/CID_SC2_selected/cid_sc2_selection.irp.f b/src/CID_SC2_selected/cid_sc2_selection.irp.f new file mode 100644 index 00000000..3610a0fe --- /dev/null +++ b/src/CID_SC2_selected/cid_sc2_selection.irp.f @@ -0,0 +1,73 @@ +program cisd_sc2_selected + implicit none + integer :: i,k + use bitmasks + + + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:),E_old(:) + integer :: N_st, iter,degree + 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 + davidson_threshold = 1.d-6 + + do while (maxval(abs(pt2(1:N_st))) > 1.d-4) + print*,'----' + print*,'' + call H_apply_cisd_selection(perturbation,pt2, norm_pert, H_pert_diag, N_st) +! soft_touch det_connections + call diagonalize_CI_SC2 + print *, 'N_det = ', N_det + do i = 1, N_st + print*,'state ',i + print *, 'PT2(SC2) = ', pt2(i) + print *, 'E(SC2) = ', CI_SC2_energy(i) + print *, 'E_before(SC2)+PT2(SC2) = ', (E_old(i)+pt2(i)) + if(i==1)then + print *, 'E(SC2)+PT2(projctd)SC2 = ', (E_old(i)+H_pert_diag(i)) + endif + E_old(i) = CI_SC2_energy(i) + enddo +! print *, 'E corr = ', (E_old(1)) - HF_energy + if (abort_all) then + exit + endif + enddo + pt2 = 0.d0 + call H_apply_PT2(pt2, norm_pert, H_pert_diag, N_st) + davidson_threshold = 1.d-10 + touch davidson_threshold davidson_criterion + do i = 1, N_st + max = 0.d0 + + print*,'' + print*,'-------------' + print*,'for state ',i + print*,'' + do k = 1, N_det + if(dabs(psi_coef(k,i)).gt.max)then + max = dabs(psi_coef(k,i)) + imax = k + endif + enddo + double precision :: max + integer :: imax + print *, 'PT2(SC2) = ', pt2(i) + print *, 'E(SC2) = ', CI_SC2_energy(i) + print *, 'E_before(SC2)+PT2(SC2) = ', (CI_SC2_energy(i)+pt2(i)) + if(i==1)then + print *, 'E(SC2)+PT2(projctd)SC2 = ', (CI_SC2_energy(i)+H_pert_diag(i)) + endif + + print*,'greater coeficient of the state : ',dabs(psi_coef(imax,i)) + call get_excitation_degree(ref_bitmask,psi_det(1,1,imax),degree,N_int) + print*,'degree of excitation of such determinant : ',degree + + enddo + print*,'coucou' + deallocate(pt2,norm_pert,H_pert_diag) +end diff --git a/src/CID_SC2_selected/do_mono_double.irp.f b/src/CID_SC2_selected/do_mono_double.irp.f new file mode 100644 index 00000000..f211879e --- /dev/null +++ b/src/CID_SC2_selected/do_mono_double.irp.f @@ -0,0 +1,19 @@ +BEGIN_PROVIDER [logical, do_double_excitations] + implicit none + BEGIN_DOC + ! if True then the double excitations are performed in the calculation + ! always true in the CISD + END_DOC + do_double_excitations = .True. + +END_PROVIDER + +BEGIN_PROVIDER [logical, do_mono_excitations] + implicit none + BEGIN_DOC + ! if True then the mono excitations are performed in the calculation + ! always true in the CISD + END_DOC + do_mono_excitations = .False. + +END_PROVIDER diff --git a/src/CID_selected/ASSUMPTIONS.rst b/src/CID_selected/ASSUMPTIONS.rst new file mode 100644 index 00000000..e69de29b diff --git a/src/CID_selected/H_apply.irp.f b/src/CID_selected/H_apply.irp.f new file mode 100644 index 00000000..91dfb9fc --- /dev/null +++ b/src/CID_selected/H_apply.irp.f @@ -0,0 +1,31 @@ +use bitmasks +BEGIN_SHELL [ /usr/bin/env python ] +from generate_h_apply import * +from perturbation import perturbations + +for perturbation in perturbations: + s = H_apply("cisd_selection_"+perturbation) + s.set_selection_pt2(perturbation) + print s +END_SHELL + + +subroutine H_apply_cisd_selection(perturbation,pt2, norm_pert, H_pert_diag, N_st) + implicit none + character*(64), intent(in) :: perturbation + integer, intent(in) :: N_st + double precision, intent(inout):: pt2(N_st), norm_pert(N_st), H_pert_diag(N_st) + +BEGIN_SHELL [ /usr/bin/env python ] +from perturbation import perturbations + +for perturbation in perturbations: + print """ + if (perturbation == '%s') then + call H_apply_cisd_selection_%s(pt2, norm_pert, H_pert_diag, N_st) + endif + """%(perturbation,perturbation) +END_SHELL + + +end diff --git a/src/CID_selected/Makefile b/src/CID_selected/Makefile new file mode 100644 index 00000000..b2ea1de1 --- /dev/null +++ b/src/CID_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/CID_selected/NEEDED_MODULES b/src/CID_selected/NEEDED_MODULES new file mode 100644 index 00000000..ac42ef59 --- /dev/null +++ b/src/CID_selected/NEEDED_MODULES @@ -0,0 +1 @@ +AOs BiInts Bitmask CISD Dets Electrons Ezfio_files Hartree_Fock MonoInts MOs Nuclei Output Perturbation SingleRefMethod Utils Selectors_full diff --git a/src/CID_selected/README.rst b/src/CID_selected/README.rst new file mode 100644 index 00000000..4fed26e5 --- /dev/null +++ b/src/CID_selected/README.rst @@ -0,0 +1,41 @@ +==================== +CISD_selected Module +==================== + +Documentation +============= + +.. Do not edit this section. It was auto-generated from the +.. NEEDED_MODULES file. + +`h_apply_cisd_selection `_ + Undocumented + +`cisd `_ + Undocumented + + + +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 `_ +* `Perturbation `_ +* `SingleRefMethod `_ +* `Utils `_ +* `Selectors_full `_ + diff --git a/src/CID_selected/cid_selection.irp.f b/src/CID_selected/cid_selection.irp.f new file mode 100644 index 00000000..f63a2d84 --- /dev/null +++ b/src/CID_selected/cid_selection.irp.f @@ -0,0 +1,34 @@ +program cisd + 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" + 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 + do i = 1, N_st + print*,'state ',i + print *, 'PT2 = ', pt2(i) + print *, 'E = ', CI_energy(i) + print *, 'E_before +PT2 = ', (E_old(i)+pt2(i)) +! print *, 'E+PT2_new= ', (E_old(1)+1.d0*pt2(1)+H_pert_diag(1))/(1.d0 +norm_pert(1)) + enddo + E_old = CI_energy + if (abort_all) then + exit + endif + enddo + deallocate(pt2,norm_pert,H_pert_diag) +end diff --git a/src/CID_selected/do_mono_double.irp.f b/src/CID_selected/do_mono_double.irp.f new file mode 100644 index 00000000..f211879e --- /dev/null +++ b/src/CID_selected/do_mono_double.irp.f @@ -0,0 +1,19 @@ +BEGIN_PROVIDER [logical, do_double_excitations] + implicit none + BEGIN_DOC + ! if True then the double excitations are performed in the calculation + ! always true in the CISD + END_DOC + do_double_excitations = .True. + +END_PROVIDER + +BEGIN_PROVIDER [logical, do_mono_excitations] + implicit none + BEGIN_DOC + ! if True then the mono excitations are performed in the calculation + ! always true in the CISD + END_DOC + do_mono_excitations = .False. + +END_PROVIDER diff --git a/src/CISD/do_mono_double.irp.f b/src/CISD/do_mono_double.irp.f new file mode 100644 index 00000000..8a07c292 --- /dev/null +++ b/src/CISD/do_mono_double.irp.f @@ -0,0 +1,19 @@ +BEGIN_PROVIDER [logical, do_double_excitations] + implicit none + BEGIN_DOC + ! if True then the double excitations are performed in the calculation + ! always true in the CISD + END_DOC + do_double_excitations = .True. + +END_PROVIDER + +BEGIN_PROVIDER [logical, do_mono_excitations] + implicit none + BEGIN_DOC + ! if True then the mono excitations are performed in the calculation + ! always true in the CISD + END_DOC + do_mono_excitations = .True. + +END_PROVIDER diff --git a/src/CISD_SC2_selected/do_mono_double.irp.f b/src/CISD_SC2_selected/do_mono_double.irp.f new file mode 100644 index 00000000..8a07c292 --- /dev/null +++ b/src/CISD_SC2_selected/do_mono_double.irp.f @@ -0,0 +1,19 @@ +BEGIN_PROVIDER [logical, do_double_excitations] + implicit none + BEGIN_DOC + ! if True then the double excitations are performed in the calculation + ! always true in the CISD + END_DOC + do_double_excitations = .True. + +END_PROVIDER + +BEGIN_PROVIDER [logical, do_mono_excitations] + implicit none + BEGIN_DOC + ! if True then the mono excitations are performed in the calculation + ! always true in the CISD + END_DOC + do_mono_excitations = .True. + +END_PROVIDER diff --git a/src/CISD_selected/do_mono_double.irp.f b/src/CISD_selected/do_mono_double.irp.f new file mode 100644 index 00000000..8a07c292 --- /dev/null +++ b/src/CISD_selected/do_mono_double.irp.f @@ -0,0 +1,19 @@ +BEGIN_PROVIDER [logical, do_double_excitations] + implicit none + BEGIN_DOC + ! if True then the double excitations are performed in the calculation + ! always true in the CISD + END_DOC + do_double_excitations = .True. + +END_PROVIDER + +BEGIN_PROVIDER [logical, do_mono_excitations] + implicit none + BEGIN_DOC + ! if True then the mono excitations are performed in the calculation + ! always true in the CISD + END_DOC + do_mono_excitations = .True. + +END_PROVIDER diff --git a/src/Dets/H_apply_template.f b/src/Dets/H_apply_template.f index 58dccb87..773fbd92 100644 --- a/src/Dets/H_apply_template.f +++ b/src/Dets/H_apply_template.f @@ -420,13 +420,17 @@ subroutine $subroutine($params_main) enddo enddo + if(do_double_excitations)then call $subroutine_diexc(psi_generators(1,1,i_generator), & mask(1,1,d_hole1), mask(1,1,d_part1), & mask(1,1,d_hole2), mask(1,1,d_part2), & i_generator $params_post) + endif + if(do_mono_excitations)then call $subroutine_monoexc(psi_generators(1,1,i_generator), & mask(1,1,s_hole ), mask(1,1,s_part ), & i_generator $params_post) + endif !$ call omp_set_lock(lck) call wall_time(wall_2) $printout_always @@ -471,13 +475,17 @@ subroutine $subroutine($params_main) not(psi_generators(k,ispin,i_generator)) ) enddo enddo - call $subroutine_diexc(psi_generators(1,1,i_generator), & - mask(1,1,d_hole1), mask(1,1,d_part1), & - mask(1,1,d_hole2), mask(1,1,d_part2), & - i_generator $params_post) - call $subroutine_monoexc(psi_generators(1,1,i_generator), & - mask(1,1,s_hole ), mask(1,1,s_part ), & - i_generator $params_post) + if(do_double_excitations)then + call $subroutine_diexc(psi_generators(1,1,i_generator), & + mask(1,1,d_hole1), mask(1,1,d_part1), & + mask(1,1,d_hole2), mask(1,1,d_part2), & + i_generator $params_post) + endif + if(do_mono_excitations)then + call $subroutine_monoexc(psi_generators(1,1,i_generator), & + mask(1,1,s_hole ), mask(1,1,s_part ), & + i_generator $params_post) + endif call wall_time(wall_2) $printout_always if (wall_2 - wall_0 > 2.d0) then diff --git a/src/Dets/filter_connected.irp.f b/src/Dets/filter_connected.irp.f index 19a1c050..e3436457 100644 --- a/src/Dets/filter_connected.irp.f +++ b/src/Dets/filter_connected.irp.f @@ -415,14 +415,14 @@ subroutine filter_connected_i_H_psi0_SC2(key1,key2,Nint,sze,idx,idx_repeat) popcnt(xor( key1(2,2,i), key2(2,2))) + & popcnt(xor( key1(3,1,i), key2(3,1))) + & popcnt(xor( key1(3,2,i), key2(3,2))) - if (degree_x2 < 5) then + if(degree_x2>6)then + idx_repeat(l_repeat) = i + l_repeat = l_repeat + 1 + else if (degree_x2 < 5) then if(degree_x2 .ne. 0)then idx(l) = i l = l+1 endif - elseif(degree_x2>6)then - idx_repeat(l_repeat) = i - l_repeat = l_repeat + 1 endif enddo From 89a7e3a6447bde2ef935e903b23ec481aa34452b Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 19 Jun 2014 17:58:45 +0200 Subject: [PATCH 3/4] DIIS on the way... --- src/Dets/README.rst | 142 ++++++++++---------- src/Hartree_Fock/Fock_matrix.irp.f | 43 ++++++ src/Hartree_Fock/HF_density_matrix_ao.irp.f | 98 ++++---------- src/Hartree_Fock/README.rst | 49 +++---- src/Hartree_Fock/SCF.irp.f | 95 +++++++++++++ src/Hartree_Fock/diagonalize_fock.irp.f | 54 ++++++-- src/Hartree_Fock/mo_SCF_iterations.irp.f | 4 +- src/MOGuess/README.rst | 8 +- src/MOs/cholesky_mo.irp.f | 80 +++++++++++ src/MOs/mos.ezfio_config | 2 + src/MOs/mos.irp.f | 22 +++ src/Perturbation/README.rst | 42 ++++-- 12 files changed, 450 insertions(+), 189 deletions(-) create mode 100644 src/Hartree_Fock/SCF.irp.f create mode 100644 src/MOs/cholesky_mo.irp.f diff --git a/src/Dets/README.rst b/src/Dets/README.rst index 35c61523..09352c44 100644 --- a/src/Dets/README.rst +++ b/src/Dets/README.rst @@ -50,24 +50,24 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`copy_h_apply_buffer_to_wf `_ +`copy_h_apply_buffer_to_wf `_ 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 `_ +`h_apply_buffer_allocated `_ Buffer of determinants/coefficients/perturbative energy for H_apply. Uninitialized. Filled by H_apply subroutines. -`h_apply_threshold `_ +`h_apply_threshold `_ Theshold on | | -`resize_h_apply_buffer `_ +`resize_h_apply_buffer `_ Undocumented -`cisd_sc2 `_ +`cisd_sc2 `_ CISD+SC2 method :: take off all the disconnected terms of a CISD (selected or not) .br dets_in : bitmasks corresponding to determinants @@ -83,29 +83,29 @@ Documentation .br Initial guess vectors are not necessarily orthonormal -`repeat_excitation `_ +`repeat_excitation `_ Undocumented -`connected_to_ref `_ +`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 -`is_in_wavefunction `_ +`is_in_wavefunction `_ Undocumented -`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 `_ +`davidson_criterion `_ Can be : [ energy | residual | both | wall_time | cpu_time | iterations ] -`davidson_diag `_ +`davidson_diag `_ Davidson diagonalization. .br dets_in : bitmasks corresponding to determinants @@ -123,7 +123,7 @@ Documentation .br Initial guess vectors are not necessarily orthonormal -`davidson_diag_hjj `_ +`davidson_diag_hjj `_ Davidson diagonalization with specific diagonal elements of the H matrix .br H_jj : specific diagonal H matrix elements to diagonalize de Davidson @@ -143,114 +143,114 @@ Documentation .br Initial guess vectors are not necessarily orthonormal -`davidson_iter_max `_ +`davidson_iter_max `_ Max number of Davidson iterations -`davidson_sze_max `_ +`davidson_sze_max `_ Max number of Davidson sizes -`davidson_threshold `_ +`davidson_threshold `_ Can be : [ energy | residual | both | wall_time | cpu_time | iterations ] -`det_search_key `_ +`det_search_key `_ Return an integer*8 corresponding to a determinant index for searching -`n_det `_ +`n_det `_ Number of determinants in the wave function -`n_det_max_jacobi `_ +`n_det_max_jacobi `_ Maximum number of determinants diagonalized my jacobi -`n_states `_ +`n_states `_ Number of states to consider -`psi_average_norm_contrib `_ +`psi_average_norm_contrib `_ Contribution of determinants to the state-averaged density -`psi_average_norm_contrib_sorted `_ +`psi_average_norm_contrib_sorted `_ Wave function sorted by determinants contribution to the norm (state-averaged) -`psi_coef `_ +`psi_coef `_ The wave function coefficients. Initialized with Hartree-Fock if the EZFIO file is empty -`psi_coef_sorted `_ +`psi_coef_sorted `_ Wave function sorted by determinants contribution to the norm (state-averaged) -`psi_coef_sorted_bit `_ +`psi_coef_sorted_bit `_ Determinants on which we apply for perturbation. o They are sorted by determinants interpreted as integers. Useful to accelerate the search of a determinant -`psi_det `_ +`psi_det `_ The wave function determinants. Initialized with Hartree-Fock if the EZFIO file is empty -`psi_det_size `_ +`psi_det_size `_ Size of the psi_det/psi_coef arrays -`psi_det_sorted `_ +`psi_det_sorted `_ Wave function sorted by determinants contribution to the norm (state-averaged) -`psi_det_sorted_bit `_ +`psi_det_sorted_bit `_ Determinants on which we apply for perturbation. o They are sorted by determinants interpreted as integers. Useful to accelerate the search of a determinant -`read_dets `_ +`read_dets `_ Reads the determinants from the EZFIO file -`save_wavefunction `_ +`save_wavefunction `_ Save the wave function into the EZFIO file -`double_exc_bitmask `_ +`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 double_exc_bitmask(:,3,i) is the bitmask for holes of excitation 2 double_exc_bitmask(:,4,i) is the bitmask for particles of excitation 2 for a given couple of hole/particle excitations i. -`n_double_exc_bitmasks `_ +`n_double_exc_bitmasks `_ Number of double excitation bitmasks -`n_single_exc_bitmasks `_ +`n_single_exc_bitmasks `_ Number of single excitation bitmasks -`single_exc_bitmask `_ +`single_exc_bitmask `_ single_exc_bitmask(:,1,i) is the bitmask for holes single_exc_bitmask(:,2,i) is the bitmask for particles for a given couple of hole/particle excitations i. -`ci_eigenvectors `_ +`ci_eigenvectors `_ Eigenvectors/values of the CI matrix -`ci_electronic_energy `_ +`ci_electronic_energy `_ Eigenvectors/values of the CI matrix -`ci_energy `_ +`ci_energy `_ N_states lowest eigenvalues of the CI matrix -`diag_algorithm `_ +`diag_algorithm `_ Diagonalization algorithm (Davidson or Lapack) -`diagonalize_ci `_ +`diagonalize_ci `_ Replace the coefficients of the CI states by the coefficients of the eigenstates of the CI matrix -`ci_sc2_eigenvectors `_ +`ci_sc2_eigenvectors `_ Eigenvectors/values of the CI matrix -`ci_sc2_electronic_energy `_ +`ci_sc2_electronic_energy `_ Eigenvectors/values of the CI matrix -`ci_sc2_energy `_ +`ci_sc2_energy `_ N_states lowest eigenvalues of the CI matrix -`diagonalize_ci_sc2 `_ +`diagonalize_ci_sc2 `_ Replace the coefficients of the CI states by the coefficients of the eigenstates of the CI matrix -`filter_connected `_ +`filter_connected `_ Filters out the determinants that are not connected by H .br returns the array idx which contains the index of the @@ -261,7 +261,7 @@ Documentation .br idx(0) is the number of determinants that interact with key1 -`filter_connected_davidson `_ +`filter_connected_davidson `_ Filters out the determinants that are not connected by H .br returns the array idx which contains the index of the @@ -272,7 +272,7 @@ Documentation .br idx(0) is the number of determinants that interact with key1 -`filter_connected_i_h_psi0 `_ +`filter_connected_i_h_psi0 `_ returns the array idx which contains the index of the .br determinants in the array key1 that interact @@ -281,7 +281,7 @@ Documentation .br idx(0) is the number of determinants that interact with key1 -`filter_connected_i_h_psi0_sc2 `_ +`filter_connected_i_h_psi0_sc2 `_ standard filter_connected_i_H_psi but returns in addition .br the array of the index of the non connected determinants to key1 @@ -292,69 +292,69 @@ Documentation .br to repeat the excitations -`get_s2 `_ +`get_s2 `_ Returns -`get_s2_u0 `_ +`get_s2_u0 `_ Undocumented -`s_z `_ +`s_z `_ Undocumented -`s_z2_sz `_ +`s_z2_sz `_ Undocumented -`a_operator `_ +`a_operator `_ Needed for diag_H_mat_elem -`ac_operator `_ +`ac_operator `_ Needed for diag_H_mat_elem -`decode_exc `_ +`decode_exc `_ Decodes the exc arrays returned by get_excitation. h1,h2 : Holes p1,p2 : Particles s1,s2 : Spins (1:alpha, 2:beta) degree : Degree of excitation -`det_connections `_ +`det_connections `_ .br -`diag_h_mat_elem `_ +`diag_h_mat_elem `_ Computes -`get_double_excitation `_ +`get_double_excitation `_ Returns the two excitation operators between two doubly excited determinants and the phase -`get_excitation `_ +`get_excitation `_ Returns the excitation operators between two determinants and the phase -`get_excitation_degree `_ +`get_excitation_degree `_ Returns the excitation degree between two determinants -`get_excitation_degree_vector `_ +`get_excitation_degree_vector `_ Applies get_excitation_degree to an array of determinants -`get_mono_excitation `_ +`get_mono_excitation `_ Returns the excitation operator between two singly excited determinants and the phase -`get_occ_from_key `_ +`get_occ_from_key `_ Returns a list of occupation numbers from a bitstring -`h_u_0 `_ +`h_u_0 `_ Computes v_0 = H|u_0> .br n : number of determinants .br H_jj : array of -`i_h_j `_ +`i_h_j `_ Returns where i and j are determinants -`i_h_psi `_ +`i_h_psi `_ for the various Nstates -`i_h_psi_sc2 `_ +`i_h_psi_sc2 `_ for the various Nstate .br returns in addition @@ -367,10 +367,10 @@ Documentation .br to repeat the excitations -`n_con_int `_ +`n_con_int `_ Number of integers to represent the connections between determinants -`h_matrix_all_dets `_ +`h_matrix_all_dets `_ H matrix on the basis of the slater deter;inants defined by psi_det diff --git a/src/Hartree_Fock/Fock_matrix.irp.f b/src/Hartree_Fock/Fock_matrix.irp.f index 42bed645..69aa7e26 100644 --- a/src/Hartree_Fock/Fock_matrix.irp.f +++ b/src/Hartree_Fock/Fock_matrix.irp.f @@ -217,3 +217,46 @@ BEGIN_PROVIDER [ double precision, HF_energy ] END_PROVIDER +BEGIN_PROVIDER [ double precision, Fock_matrix_ao, (ao_num_align, ao_num) ] + implicit none + BEGIN_DOC + ! Fock matrix in AO basis set + END_DOC + + if (elec_alpha_num == elec_beta_num) then + integer :: i,j + do j=1,ao_num + !DIR$ VECTOR ALIGNED + do i=1,ao_num_align + Fock_matrix_ao(i,j) = Fock_matrix_alpha_ao(i,j) + enddo + enddo + else + double precision, allocatable :: T(:,:), M(:,:) + ! F_ao = S C F_mo C^t S + allocate (T(mo_tot_num_align,mo_tot_num),M(ao_num_align,mo_tot_num)) + call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & + ao_overlap, size(ao_overlap,1), & + mo_coef, size(mo_coef,1), & + 0.d0, & + M, size(M,1)) + call dgemm('N','N', ao_num,mo_tot_num,mo_tot_num, 1.d0, & + M, size(M,1), & + Fock_matrix_mo, size(Fock_matrix_mo,1), & + 0.d0, & + T, size(T,1)) + call dgemm('N','T', mo_tot_num,ao_num,mo_tot_num, 1.d0, & + T, size(T,1), & + mo_coef, size(mo_coef,1), & + 0.d0, & + M, size(M,1)) + call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & + M, size(M,1), & + ao_overlap, size(ao_overlap,1), & + 0.d0, & + Fock_matrix_ao, size(Fock_matrix_ao,1)) + + deallocate(T) + endif +END_PROVIDER + diff --git a/src/Hartree_Fock/HF_density_matrix_ao.irp.f b/src/Hartree_Fock/HF_density_matrix_ao.irp.f index b0d58344..e85e4a6c 100644 --- a/src/Hartree_Fock/HF_density_matrix_ao.irp.f +++ b/src/Hartree_Fock/HF_density_matrix_ao.irp.f @@ -1,46 +1,27 @@ - BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_alpha, (ao_num_align,ao_num) ] -&BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_beta, (ao_num_align,ao_num) ] +BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_alpha, (ao_num_align,ao_num) ] implicit none BEGIN_DOC - ! Alpha and Beta density matrix in the AO basis + ! Alpha density matrix in the AO basis END_DOC - integer :: i,j,k,l1,l2 - integer, allocatable :: mo_occ(:,:) - allocate ( mo_occ(elec_alpha_num,2) ) - call bitstring_to_list( HF_bitmask(1,1), mo_occ(1,1), j, N_int) - ASSERT ( j==elec_alpha_num ) + call dgemm('N','T',ao_num,ao_num,elec_alpha_num,1.d0, & + mo_coef, size(mo_coef,1), & + mo_coef, size(mo_coef,1), 0.d0, & + HF_density_matrix_ao_alpha, size(HF_density_matrix_ao_alpha,1)) + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_beta, (ao_num_align,ao_num) ] + implicit none + BEGIN_DOC + ! Beta density matrix in the AO basis + END_DOC - call bitstring_to_list( HF_bitmask(1,2), mo_occ(1,2), j, N_int) - ASSERT ( j==elec_beta_num ) - - do j=1,ao_num - !DIR$ VECTOR ALIGNED - do i=1,ao_num_align - HF_density_matrix_ao_alpha(i,j) = 0.d0 - HF_density_matrix_ao_beta (i,j) = 0.d0 - enddo - do k=1,elec_beta_num - l1 = mo_occ(k,1) - l2 = mo_occ(k,2) - !DIR$ VECTOR ALIGNED - do i=1,ao_num - HF_density_matrix_ao_alpha(i,j) = HF_density_matrix_ao_alpha(i,j) +& - mo_coef(i,l1) * mo_coef(j,l1) - HF_density_matrix_ao_beta (i,j) = HF_density_matrix_ao_beta (i,j) +& - mo_coef(i,l2) * mo_coef(j,l2) - enddo - enddo - do k=elec_beta_num+1,elec_alpha_num - l1 = mo_occ(k,1) - !DIR$ VECTOR ALIGNED - do i=1,ao_num - HF_density_matrix_ao_alpha(i,j) = HF_density_matrix_ao_alpha(i,j) +& - mo_coef(i,l1) * mo_coef(j,l1) - enddo - enddo - enddo - deallocate(mo_occ) + call dgemm('N','T',ao_num,ao_num,elec_beta_num,1.d0, & + mo_coef, size(mo_coef,1), & + mo_coef, size(mo_coef,1), 0.d0, & + HF_density_matrix_ao_beta, size(HF_density_matrix_ao_beta,1)) + END_PROVIDER BEGIN_PROVIDER [ double precision, HF_density_matrix_ao, (ao_num_align,ao_num) ] @@ -48,40 +29,13 @@ BEGIN_PROVIDER [ double precision, HF_density_matrix_ao, (ao_num_align,ao_num) ] BEGIN_DOC ! Density matrix in the AO basis END_DOC - integer :: i,j,k,l1,l2 - integer, allocatable :: mo_occ(:,:) + ASSERT (size(HF_density_matrix_ao,1) == size(HF_density_matrix_ao_alpha,1)) + if (elec_alpha_num== elec_beta_num) then + HF_density_matrix_ao = HF_density_matrix_ao_alpha + HF_density_matrix_ao_alpha + else + ASSERT (size(HF_density_matrix_ao,1) == size(HF_density_matrix_ao_beta ,1)) + HF_density_matrix_ao = HF_density_matrix_ao_alpha + HF_density_matrix_ao_beta + endif - allocate ( mo_occ(elec_alpha_num,2) ) - call bitstring_to_list( HF_bitmask(1,1), mo_occ(1,1), j, N_int) - ASSERT ( j==elec_alpha_num ) - - call bitstring_to_list( HF_bitmask(1,2), mo_occ(1,2), j, N_int) - ASSERT ( j==elec_beta_num ) - - do j=1,ao_num - !DIR$ VECTOR ALIGNED - do i=1,ao_num_align - HF_density_matrix_ao(i,j) = 0.d0 - enddo - do k=1,elec_beta_num - l1 = mo_occ(k,1) - l2 = mo_occ(k,2) - !DIR$ VECTOR ALIGNED - do i=1,ao_num - HF_density_matrix_ao(i,j) = HF_density_matrix_ao(i,j) + & - mo_coef(i,l1) * mo_coef(j,l1) + & - mo_coef(i,l2) * mo_coef(j,l2) - enddo - enddo - do k=elec_beta_num+1,elec_alpha_num - l1 = mo_occ(k,1) - !DIR$ VECTOR ALIGNED - do i=1,ao_num - HF_density_matrix_ao(i,j) = HF_density_matrix_ao(i,j) + & - mo_coef(i,l1) * mo_coef(j,l1) - enddo - enddo - enddo - deallocate(mo_occ) END_PROVIDER diff --git a/src/Hartree_Fock/README.rst b/src/Hartree_Fock/README.rst index 074928ad..7173d418 100644 --- a/src/Hartree_Fock/README.rst +++ b/src/Hartree_Fock/README.rst @@ -26,19 +26,22 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`fock_matrix_alpha_ao `_ +`fock_matrix_alpha_ao `_ Alpha Fock matrix in AO basis set -`fock_matrix_alpha_mo `_ +`fock_matrix_alpha_mo `_ Fock matrix on the MO basis -`fock_matrix_beta_ao `_ +`fock_matrix_ao `_ + Fock matrix in AO basis set + +`fock_matrix_beta_ao `_ Alpha Fock matrix in AO basis set -`fock_matrix_beta_mo `_ +`fock_matrix_beta_mo `_ Fock matrix on the MO basis -`fock_matrix_diag_mo `_ +`fock_matrix_diag_mo `_ Fock matrix on the MO basis. For open shells, the ROHF Fock Matrix is .br @@ -53,7 +56,7 @@ Documentation K = Fb - Fa .br -`fock_matrix_mo `_ +`fock_matrix_mo `_ Fock matrix on the MO basis. For open shells, the ROHF Fock Matrix is .br @@ -68,49 +71,49 @@ Documentation K = Fb - Fa .br -`hf_energy `_ +`hf_energy `_ Hartree-Fock energy -`hf_density_matrix_ao `_ +`hf_density_matrix_ao `_ Density matrix in the AO basis -`hf_density_matrix_ao_alpha `_ - Alpha and Beta density matrix in the AO basis +`hf_density_matrix_ao_alpha `_ + Alpha density matrix in the AO basis -`hf_density_matrix_ao_beta `_ - Alpha and Beta density matrix in the AO basis +`hf_density_matrix_ao_beta `_ + Beta density matrix in the AO basis -`diagonal_fock_matrix_mo `_ +`diagonal_fock_matrix_mo `_ Diagonal Fock matrix in the MO basis -`eigenvectors_fock_matrix_mo `_ +`eigenvectors_fock_matrix_mo `_ Diagonal Fock matrix in the MO basis -`scf_iteration `_ +`xcf_iteration `_ Undocumented -`do_diis `_ +`do_diis `_ If True, compute integrals on the fly -`n_it_scf_max `_ +`n_it_scf_max `_ Maximum number of SCF iterations -`thresh_scf `_ +`thresh_scf `_ Threshold on the convergence of the Hartree Fock energy -`bi_elec_ref_bitmask_energy `_ +`bi_elec_ref_bitmask_energy `_ Energy of the reference bitmask used in Slater rules -`kinetic_ref_bitmask_energy `_ +`kinetic_ref_bitmask_energy `_ Energy of the reference bitmask used in Slater rules -`mono_elec_ref_bitmask_energy `_ +`mono_elec_ref_bitmask_energy `_ Energy of the reference bitmask used in Slater rules -`nucl_elec_ref_bitmask_energy `_ +`nucl_elec_ref_bitmask_energy `_ Energy of the reference bitmask used in Slater rules -`ref_bitmask_energy `_ +`ref_bitmask_energy `_ Energy of the reference bitmask used in Slater rules diff --git a/src/Hartree_Fock/SCF.irp.f b/src/Hartree_Fock/SCF.irp.f new file mode 100644 index 00000000..cb6e8cbf --- /dev/null +++ b/src/Hartree_Fock/SCF.irp.f @@ -0,0 +1,95 @@ +BEGIN_PROVIDER [ integer, it_scf ] + implicit none + BEGIN_DOC + ! Number of the current SCF iteration + END_DOC + it_scf = 0 +END_PROVIDER + +BEGIN_PROVIDER [ double precision, SCF_density_matrices, (ao_num_align,ao_num,2,0:n_it_scf_max) ] + implicit none + BEGIN_DOC + ! Density matrices at every SCF iteration + END_DOC + SCF_density_matrices = 0.d0 +END_PROVIDER + +subroutine insert_new_SCF_density_matrix + implicit none + integer :: i,j + do j=1,ao_num + do i=1,ao_num + SCF_density_matrices(i,j,1,it_scf) = HF_density_matrix_ao_alpha(i,j) + SCF_density_matrices(i,j,1,0) += HF_density_matrix_ao_alpha(i,j) + SCF_density_matrices(i,j,2,it_scf) = HF_density_matrix_ao_beta(i,j) + SCF_density_matrices(i,j,2,0) += HF_density_matrix_ao_beta(i,j) + enddo + enddo +end + +subroutine Fock_mo_to_ao(FMO,LDFMO,FAO,LDFAO) + implicit none + integer, intent(in) :: LDFMO ! size(FMO,1) + integer, intent(in) :: LDFAO ! size(FAO,1) + double precision, intent(in) :: FMO(LDFMO,*) + double precision, intent(out) :: FAO(LDFAO,*) + + double precision, allocatable :: T(:,:), M(:,:) + ! F_ao = S C F_mo C^t S + allocate (T(mo_tot_num_align,mo_tot_num),M(ao_num_align,mo_tot_num)) + call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & + ao_overlap, size(ao_overlap,1), & + mo_coef, size(mo_coef,1), & + 0.d0, & + M, size(M,1)) + call dgemm('N','N', ao_num,mo_tot_num,mo_tot_num, 1.d0, & + M, size(M,1), & + FMO, size(FMO,1), & + 0.d0, & + T, size(T,1)) + call dgemm('N','T', mo_tot_num,ao_num,mo_tot_num, 1.d0, & + T, size(T,1), & + mo_coef, size(mo_coef,1), & + 0.d0, & + M, size(M,1)) + call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & + M, size(M,1), & + ao_overlap, size(ao_overlap,1), & + 0.d0, & + FAO, size(FAO,1)) + deallocate(T,M) +end + +subroutine DIIS_step + implicit none + integer :: i,j + double precision :: c + c = 1.d0/dble(it_scf) + do j=1,ao_num + do i=1,ao_num + HF_density_matrix_ao_alpha(i,j) = SCF_density_matrices(i,j,1,0) * c + HF_density_matrix_ao_beta (i,j) = SCF_density_matrices(i,j,2,0) * c + enddo + enddo + TOUCH HF_density_matrix_ao_alpha HF_density_matrix_ao_beta + +! call Fock_mo_to_ao(Fock_matrix_mo_alpha, size(Fock_matrix_mo_alpha,1), & +! Fock_matrix_alpha_ao, size(Fock_matrix_alpha_ao,1) ) +! call Fock_mo_to_ao(Fock_matrix_mo_beta, size(Fock_matrix_mo_beta,1), & +! Fock_matrix_beta_ao, size(Fock_matrix_beta_ao,1) ) +! SOFT_TOUCH Fock_matrix_alpha_ao Fock_matrix_beta_ao Fock_matrix_mo_alpha Fock_matrix_mo_beta +end + +subroutine scf_iteration + implicit none + integer :: i,j + do i=1,n_it_scf_max + it_scf += 1 + SOFT_TOUCH it_scf + mo_coef = eigenvectors_Fock_matrix_mo + TOUCH mo_coef + call insert_new_SCF_density_matrix + call DIIS_step + print *, HF_energy + enddo +end diff --git a/src/Hartree_Fock/diagonalize_fock.irp.f b/src/Hartree_Fock/diagonalize_fock.irp.f index c34a1415..c5faea9c 100644 --- a/src/Hartree_Fock/diagonalize_fock.irp.f +++ b/src/Hartree_Fock/diagonalize_fock.irp.f @@ -5,16 +5,52 @@ ! Diagonal Fock matrix in the MO basis END_DOC - double precision, allocatable :: R(:,:) - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: R + integer :: i,j + integer :: liwork, lwork, n, info + integer, allocatable :: iwork(:) + double precision, allocatable :: work(:), F(:,:), S(:,:) - allocate(R(mo_tot_num_align,mo_tot_num)) + allocate(F(ao_num_align,ao_num), S(ao_num_align,ao_num) ) + do j=1,ao_num + do i=1,ao_num + S(i,j) = ao_overlap(i,j) + F(i,j) = Fock_matrix_ao(i,j) + enddo + enddo + + n = ao_num + lwork = 1+6*n + 2*n*n + liwork = 3 + 5*n - call lapack_diag(diagonal_Fock_matrix_mo,R,Fock_matrix_mo,size(Fock_matrix_mo,1),& - mo_tot_num) - - call dgemm('N','N',ao_num,mo_tot_num,mo_tot_num,1.d0,mo_coef,size(mo_coef,1),& - R,size(R,1),0.d0,eigenvectors_Fock_matrix_mo,size(eigenvectors_Fock_matrix_mo,1)) - deallocate(R) + allocate(work(lwork), iwork(liwork) ) + + lwork = -1 + liwork = -1 + call dsygvd(1,'v','u',mo_tot_num,F,size(F,1),S,size(S,1),& + diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info) + + if (info /= 0) then + print *, irp_here//' failed' + stop 1 + endif + lwork = int(work(1)) + liwork = iwork(1) + deallocate(work,iwork) + allocate(work(lwork), iwork(liwork) ) + + call dsygvd(1,'v','u',mo_tot_num,F,size(F,1),S,size(S,1),& + diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info) + + if (info /= 0) then + print *, irp_here//' failed' + stop 1 + endif + do j=1,ao_num + do i=1,ao_num + eigenvectors_Fock_matrix_mo(i,j) = F(i,j) + enddo + enddo + + deallocate(work, iwork, F, S) END_PROVIDER diff --git a/src/Hartree_Fock/mo_SCF_iterations.irp.f b/src/Hartree_Fock/mo_SCF_iterations.irp.f index 01cab413..4854861f 100644 --- a/src/Hartree_Fock/mo_SCF_iterations.irp.f +++ b/src/Hartree_Fock/mo_SCF_iterations.irp.f @@ -1,4 +1,4 @@ -program scf_iteration +program xcf_iteration use bitmasks implicit none double precision :: SCF_energy_before,SCF_energy_after,diag_H_mat_elem,get_mo_bielec_integral @@ -32,6 +32,6 @@ program scf_iteration endif mo_label = "Canonical" TOUCH mo_label mo_coef - call save_mos +! call save_mos end diff --git a/src/MOGuess/README.rst b/src/MOGuess/README.rst index d8e72641..90785358 100644 --- a/src/MOGuess/README.rst +++ b/src/MOGuess/README.rst @@ -22,19 +22,19 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`h_core_guess `_ +`h_core_guess `_ Undocumented -`ao_ortho_lowdin_coef `_ +`ao_ortho_lowdin_coef `_ matrix of the coefficients of the mos generated by the orthonormalization by the S^{-1/2} canonical transformation of the aos ao_ortho_lowdin_coef(i,j) = coefficient of the ith ao on the jth ao_ortho_lowdin orbital -`ao_ortho_lowdin_overlap `_ +`ao_ortho_lowdin_overlap `_ overlap matrix of the ao_ortho_lowdin supposed to be the Identity -`ao_ortho_lowdin_nucl_elec_integral `_ +`ao_ortho_lowdin_nucl_elec_integral `_ Undocumented diff --git a/src/MOs/cholesky_mo.irp.f b/src/MOs/cholesky_mo.irp.f new file mode 100644 index 00000000..97b6abd2 --- /dev/null +++ b/src/MOs/cholesky_mo.irp.f @@ -0,0 +1,80 @@ +subroutine cholesky_mo(n,m,P,LDP,C,LDC,tol_in,rank) + implicit none + BEGIN_DOC +! Cholesky decomposition of AO Density matrix to +! generate MOs + END_DOC + integer, intent(in) :: n,m, LDC, LDP + double precision, intent(in) :: P(LDP,n) + double precision, intent(out) :: C(LDC,m) + double precision, intent(in) :: tol_in + integer, intent(out) :: rank + + integer :: info + integer :: i,k + integer :: ipiv(n) + double precision:: tol + double precision, allocatable :: W(:,:), work(:) + !DEC$ ATTRIBUTES ALIGN: 32 :: W + !DEC$ ATTRIBUTES ALIGN: 32 :: work + !DEC$ ATTRIBUTES ALIGN: 32 :: ipiv + + allocate(W(LDC,n),work(2*n)) + tol=tol_in + + info = 0 + do i=1,n + do k=1,i + W(i,k) = P(i,k) + enddo + do k=i+1,n + W(i,k) = 0. + enddo + enddo + call DPSTRF('L', n, W, LDC, ipiv, rank, tol, work, info ) + do i=1,n + do k=1,min(m,rank) + C(ipiv(i),k) = W(i,k) + enddo + enddo + + deallocate(W,work) +end + +BEGIN_PROVIDER [ double precision, mo_density_matrix, (mo_tot_num_align, mo_tot_num) ] + implicit none + BEGIN_DOC + ! Density matrix in MO basis + END_DOC + integer :: i,j,k + mo_density_matrix = 0.d0 + do k=1,mo_tot_num + if (mo_occ(k) == 0.d0) then + cycle + endif + do j=1,ao_num + do i=1,ao_num + mo_density_matrix(i,j) = mo_density_matrix(i,j) + & + mo_occ(k) * mo_coef(i,k) * mo_coef(j,k) + enddo + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [ double precision, mo_density_matrix_virtual, (mo_tot_num_align, mo_tot_num) ] + implicit none + BEGIN_DOC + ! Density matrix in MO basis (virtual MOs) + END_DOC + integer :: i,j,k + mo_density_matrix_virtual = 0.d0 + do k=1,mo_tot_num + do j=1,ao_num + do i=1,ao_num + mo_density_matrix_virtual(i,j) = mo_density_matrix_virtual(i,j) + & + (2.d0-mo_occ(k)) * mo_coef(i,k) * mo_coef(j,k) + enddo + enddo + enddo +END_PROVIDER + diff --git a/src/MOs/mos.ezfio_config b/src/MOs/mos.ezfio_config index a0eda491..e34d1aab 100644 --- a/src/MOs/mos.ezfio_config +++ b/src/MOs/mos.ezfio_config @@ -1,5 +1,7 @@ mo_basis + ao_num integer mo_tot_num integer mo_coef double precision (ao_basis_ao_num,mo_basis_mo_tot_num) mo_label character*(64) + mo_occ double precision (mo_basis_mo_tot_num) diff --git a/src/MOs/mos.irp.f b/src/MOs/mos.irp.f index b081d8ce..7895ac55 100644 --- a/src/MOs/mos.irp.f +++ b/src/MOs/mos.irp.f @@ -75,3 +75,25 @@ BEGIN_PROVIDER [ double precision, mo_coef_transp, (mo_tot_num_align,ao_num) ] END_PROVIDER +BEGIN_PROVIDER [ double precision, mo_occ, (mo_tot_num) ] + implicit none + BEGIN_DOC + ! MO occupation numbers + END_DOC + PROVIDE ezfio_filename + logical :: exists + call ezfio_has_mo_basis_mo_occ(exists) + if (exists) then + call ezfio_get_mo_basis_mo_occ(mo_occ) + else + mo_occ = 0.d0 + integer :: i + do i=1,elec_beta_num + mo_occ(i) = 2.d0 + enddo + do i=elec_beta_num+1,elec_alpha_num + mo_occ(i) = 1.d0 + enddo + endif +END_PROVIDER + diff --git a/src/Perturbation/README.rst b/src/Perturbation/README.rst index 806ad292..aad1d0c0 100644 --- a/src/Perturbation/README.rst +++ b/src/Perturbation/README.rst @@ -82,7 +82,7 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`pt2_moller_plesset `_ +`pt2_moller_plesset `_ compute the standard Moller-Plesset perturbative first order coefficient and second order energetic contribution .br for the various n_st states. @@ -92,7 +92,7 @@ Documentation e_2_pert(i) = ^2/(difference of orbital energies) .br -`pt2_epstein_nesbet `_ +`pt2_epstein_nesbet `_ compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution .br for the various N_st states. @@ -102,7 +102,7 @@ Documentation e_2_pert(i) = ^2/( E(i) - ) .br -`pt2_epstein_nesbet_2x2 `_ +`pt2_epstein_nesbet_2x2 `_ compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution .br for the various N_st states. @@ -112,20 +112,46 @@ Documentation c_pert(i) = e_2_pert(i)/ .br -`fill_h_apply_buffer_selection `_ +`pt2_epstein_nesbet_sc2_projected `_ + compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution + .br + for the various N_st states, + .br + but with the correction in the denominator + .br + comming from the interaction of that determinant with all the others determinants + .br + that can be repeated by repeating all the double excitations + .br + : you repeat all the correlation energy already taken into account in CI_electronic_energy(1) + .br + that could be repeated to this determinant. + .br + In addition, for the perturbative energetic contribution you have the standard second order + .br + e_2_pert = ^2/(Delta_E) + .br + and also the purely projected contribution + .br + H_pert_diag = c_pert + +`repeat_all_e_corr `_ + Undocumented + +`fill_h_apply_buffer_selection `_ Fill the H_apply buffer with determiants for the selection -`remove_small_contributions `_ +`remove_small_contributions `_ Remove determinants with small contributions. N_states is assumed to be provided. -`selection_criterion `_ +`selection_criterion `_ Threshold to select determinants. Set by selection routines. -`selection_criterion_factor `_ +`selection_criterion_factor `_ Threshold to select determinants. Set by selection routines. -`selection_criterion_min `_ +`selection_criterion_min `_ Threshold to select determinants. Set by selection routines. From 12c47364ca83852d01a1e6c550c57f385d8a17d7 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 19 Jun 2014 22:34:56 +0200 Subject: [PATCH 4/4] Better Hartree-Fock --- src/Hartree_Fock/README.rst | 21 +++ src/Hartree_Fock/SCF.irp.f | 170 ++++++++++++----------- src/Hartree_Fock/mo_SCF_iterations.irp.f | 26 +--- 3 files changed, 114 insertions(+), 103 deletions(-) diff --git a/src/Hartree_Fock/README.rst b/src/Hartree_Fock/README.rst index 7173d418..a9d70bd4 100644 --- a/src/Hartree_Fock/README.rst +++ b/src/Hartree_Fock/README.rst @@ -83,6 +83,27 @@ Documentation `hf_density_matrix_ao_beta `_ Beta density matrix in the AO basis +`fock_mo_to_ao `_ + Undocumented + +`insert_new_scf_density_matrix `_ + Undocumented + +`it_scf `_ + Number of the current SCF iteration + +`scf_density_matrices `_ + Density matrices at every SCF iteration + +`scf_energies `_ + Density matrices at every SCF iteration + +`scf_interpolation_step `_ + Undocumented + +`scf_iterations `_ + Undocumented + `diagonal_fock_matrix_mo `_ Diagonal Fock matrix in the MO basis diff --git a/src/Hartree_Fock/SCF.irp.f b/src/Hartree_Fock/SCF.irp.f index cb6e8cbf..895ca546 100644 --- a/src/Hartree_Fock/SCF.irp.f +++ b/src/Hartree_Fock/SCF.irp.f @@ -1,95 +1,107 @@ BEGIN_PROVIDER [ integer, it_scf ] - implicit none - BEGIN_DOC - ! Number of the current SCF iteration - END_DOC - it_scf = 0 + implicit none + BEGIN_DOC + ! Number of the current SCF iteration + END_DOC + it_scf = 0 END_PROVIDER -BEGIN_PROVIDER [ double precision, SCF_density_matrices, (ao_num_align,ao_num,2,0:n_it_scf_max) ] - implicit none - BEGIN_DOC - ! Density matrices at every SCF iteration - END_DOC - SCF_density_matrices = 0.d0 + BEGIN_PROVIDER [ double precision, SCF_density_matrices, (ao_num_align,ao_num,2,n_it_scf_max) ] +&BEGIN_PROVIDER [ double precision, SCF_energies, (n_it_scf_max) ] + implicit none + BEGIN_DOC + ! Density matrices at every SCF iteration + END_DOC + SCF_density_matrices = 0.d0 + SCF_energies = 0.d0 END_PROVIDER subroutine insert_new_SCF_density_matrix - implicit none - integer :: i,j - do j=1,ao_num - do i=1,ao_num - SCF_density_matrices(i,j,1,it_scf) = HF_density_matrix_ao_alpha(i,j) - SCF_density_matrices(i,j,1,0) += HF_density_matrix_ao_alpha(i,j) - SCF_density_matrices(i,j,2,it_scf) = HF_density_matrix_ao_beta(i,j) - SCF_density_matrices(i,j,2,0) += HF_density_matrix_ao_beta(i,j) + implicit none + integer :: i,j + do j=1,ao_num + do i=1,ao_num + SCF_density_matrices(i,j,1,it_scf) = HF_density_matrix_ao_alpha(i,j) + SCF_density_matrices(i,j,2,it_scf) = HF_density_matrix_ao_beta(i,j) + enddo enddo - enddo + SCF_energies(it_scf) = HF_energy end subroutine Fock_mo_to_ao(FMO,LDFMO,FAO,LDFAO) - implicit none - integer, intent(in) :: LDFMO ! size(FMO,1) - integer, intent(in) :: LDFAO ! size(FAO,1) - double precision, intent(in) :: FMO(LDFMO,*) - double precision, intent(out) :: FAO(LDFAO,*) - - double precision, allocatable :: T(:,:), M(:,:) - ! F_ao = S C F_mo C^t S - allocate (T(mo_tot_num_align,mo_tot_num),M(ao_num_align,mo_tot_num)) - call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & - ao_overlap, size(ao_overlap,1), & - mo_coef, size(mo_coef,1), & - 0.d0, & - M, size(M,1)) - call dgemm('N','N', ao_num,mo_tot_num,mo_tot_num, 1.d0, & - M, size(M,1), & - FMO, size(FMO,1), & - 0.d0, & - T, size(T,1)) - call dgemm('N','T', mo_tot_num,ao_num,mo_tot_num, 1.d0, & - T, size(T,1), & - mo_coef, size(mo_coef,1), & - 0.d0, & - M, size(M,1)) - call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & - M, size(M,1), & - ao_overlap, size(ao_overlap,1), & - 0.d0, & - FAO, size(FAO,1)) - deallocate(T,M) + implicit none + integer, intent(in) :: LDFMO ! size(FMO,1) + integer, intent(in) :: LDFAO ! size(FAO,1) + double precision, intent(in) :: FMO(LDFMO,*) + double precision, intent(out) :: FAO(LDFAO,*) + + double precision, allocatable :: T(:,:), M(:,:) + ! F_ao = S C F_mo C^t S + allocate (T(mo_tot_num_align,mo_tot_num),M(ao_num_align,mo_tot_num)) + call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & + ao_overlap, size(ao_overlap,1), & + mo_coef, size(mo_coef,1), & + 0.d0, & + M, size(M,1)) + call dgemm('N','N', ao_num,mo_tot_num,mo_tot_num, 1.d0, & + M, size(M,1), & + FMO, size(FMO,1), & + 0.d0, & + T, size(T,1)) + call dgemm('N','T', mo_tot_num,ao_num,mo_tot_num, 1.d0, & + T, size(T,1), & + mo_coef, size(mo_coef,1), & + 0.d0, & + M, size(M,1)) + call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & + M, size(M,1), & + ao_overlap, size(ao_overlap,1), & + 0.d0, & + FAO, size(FAO,1)) + deallocate(T,M) end -subroutine DIIS_step - implicit none - integer :: i,j - double precision :: c - c = 1.d0/dble(it_scf) - do j=1,ao_num - do i=1,ao_num - HF_density_matrix_ao_alpha(i,j) = SCF_density_matrices(i,j,1,0) * c - HF_density_matrix_ao_beta (i,j) = SCF_density_matrices(i,j,2,0) * c +subroutine SCF_interpolation_step + implicit none + integer :: i,j + double precision :: c + + if (it_scf == 1) then + return + endif + call random_number(c) + do j=1,ao_num + do i=1,ao_num + HF_density_matrix_ao_alpha(i,j) = c*SCF_density_matrices(i,j,1,it_scf)+SCF_density_matrices(i,j,1,it_scf-1) * (1.d0 - c) + HF_density_matrix_ao_beta (i,j) = c*SCF_density_matrices(i,j,2,it_scf)+SCF_density_matrices(i,j,2,it_scf-1) * (1.d0 - c) + enddo enddo - enddo - TOUCH HF_density_matrix_ao_alpha HF_density_matrix_ao_beta - -! call Fock_mo_to_ao(Fock_matrix_mo_alpha, size(Fock_matrix_mo_alpha,1), & -! Fock_matrix_alpha_ao, size(Fock_matrix_alpha_ao,1) ) -! call Fock_mo_to_ao(Fock_matrix_mo_beta, size(Fock_matrix_mo_beta,1), & -! Fock_matrix_beta_ao, size(Fock_matrix_beta_ao,1) ) -! SOFT_TOUCH Fock_matrix_alpha_ao Fock_matrix_beta_ao Fock_matrix_mo_alpha Fock_matrix_mo_beta + TOUCH HF_density_matrix_ao_alpha HF_density_matrix_ao_beta + + ! call Fock_mo_to_ao(Fock_matrix_mo_alpha, size(Fock_matrix_mo_alpha,1),& + ! Fock_matrix_alpha_ao, size(Fock_matrix_alpha_ao,1) ) + ! call Fock_mo_to_ao(Fock_matrix_mo_beta, size(Fock_matrix_mo_beta,1),& + ! Fock_matrix_beta_ao, size(Fock_matrix_beta_ao,1) ) + ! SOFT_TOUCH Fock_matrix_alpha_ao Fock_matrix_beta_ao Fock_matrix_mo_alpha Fock_matrix_mo_beta end -subroutine scf_iteration - implicit none - integer :: i,j - do i=1,n_it_scf_max - it_scf += 1 - SOFT_TOUCH it_scf - mo_coef = eigenvectors_Fock_matrix_mo - TOUCH mo_coef - call insert_new_SCF_density_matrix - call DIIS_step - print *, HF_energy - enddo +subroutine scf_iterations + implicit none + integer :: i,j + do i=1,n_it_scf_max + it_scf += 1 + SOFT_TOUCH it_scf + mo_coef = eigenvectors_Fock_matrix_mo + TOUCH mo_coef + call insert_new_SCF_density_matrix + print *, HF_energy + if (SCF_energies(it_scf)>SCF_energies(it_scf-1)) then + call SCF_interpolation_step + endif + if (it_scf>1 ) then + if (dabs(SCF_energies(it_scf)-SCF_energies(it_scf-1)) < thresh_SCF) then + exit + endif + endif + enddo end diff --git a/src/Hartree_Fock/mo_SCF_iterations.irp.f b/src/Hartree_Fock/mo_SCF_iterations.irp.f index 4854861f..bae0e6d4 100644 --- a/src/Hartree_Fock/mo_SCF_iterations.irp.f +++ b/src/Hartree_Fock/mo_SCF_iterations.irp.f @@ -6,32 +6,10 @@ program xcf_iteration integer :: i_it E0 = HF_energy - i_it = 0 - n_it_scf_max = 10 - SCF_energy_before = huge(1.d0) - SCF_energy_after = E0 - print *, E0 - mo_label = "Canonical" thresh_SCF = 1.d-10 - do while (i_it < 40 .and. dabs(SCF_energy_before - SCF_energy_after) > thresh_SCF) - SCF_energy_before = SCF_energy_after - mo_coef = eigenvectors_Fock_matrix_mo - TOUCH mo_coef mo_label - call clear_mo_map - SCF_energy_after = HF_energy - print*,SCF_energy_after, dabs(SCF_energy_before - SCF_energy_after) - i_it +=1 - if(i_it > n_it_scf_max)exit - enddo - - if (i_it >= n_it_scf_max) then - stop 'Failed' - endif - if (SCF_energy_after - E0 > thresh_SCF) then - stop 'Failed' - endif + call scf_iterations mo_label = "Canonical" TOUCH mo_label mo_coef -! call save_mos + call save_mos end