add CID, CID_selected, CID_SC2_selected, just for fun ...

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
+

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

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
+

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 <http://github.com/LCPQ/quantum_package/tree/master/src/AOs>`_
+* `BiInts <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts>`_
+* `Bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask>`_
+* `Dets <http://github.com/LCPQ/quantum_package/tree/master/src/Dets>`_
+* `Electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Electrons>`_
+* `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
+* `Hartree_Fock <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock>`_
+* `MonoInts <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts>`_
+* `MOs <http://github.com/LCPQ/quantum_package/tree/master/src/MOs>`_
+* `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_
+* `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
+* `SingleRefMethod <http://github.com/LCPQ/quantum_package/tree/master/src/SingleRefMethod>`_
+* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
+* `Selectors_full <http://github.com/LCPQ/quantum_package/tree/master/src/Selectors_full>`_
+
+Documentation
+=============
+
+.. Do not edit this section. It was auto-generated from the
+.. NEEDED_MODULES file.
+
+`cisd <http://github.com/LCPQ/quantum_package/tree/master/src/CISD/cisd_lapack.irp.f#L1>`_
+  Undocumented
+
+
+

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

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

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

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
+

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

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

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 <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_SC2_selected/cisd_sc2_selection.irp.f#L1>`_
+  Undocumented
+
+
+
+Needed Modules
+==============
+
+.. Do not edit this section. It was auto-generated from the
+.. NEEDED_MODULES file.
+
+* `AOs <http://github.com/LCPQ/quantum_package/tree/master/src/AOs>`_
+* `BiInts <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts>`_
+* `Bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask>`_
+* `CISD <http://github.com/LCPQ/quantum_package/tree/master/src/CISD>`_
+* `SC2 <http://github.com/LCPQ/quantum_package/tree/master/src/SC2>`_
+* `CISD_selected <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_selected>`_
+* `Dets <http://github.com/LCPQ/quantum_package/tree/master/src/Dets>`_
+* `Electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Electrons>`_
+* `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
+* `Hartree_Fock <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock>`_
+* `MonoInts <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts>`_
+* `MOs <http://github.com/LCPQ/quantum_package/tree/master/src/MOs>`_
+* `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_
+* `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
+* `Perturbation <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation>`_
+* `Selectors_full <http://github.com/LCPQ/quantum_package/tree/master/src/Selectors_full>`_
+* `SingleRefMethod <http://github.com/LCPQ/quantum_package/tree/master/src/SingleRefMethod>`_
+* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
+

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

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

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

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

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

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 <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_selected/H_apply.irp.f#L13>`_
+  Undocumented
+
+`cisd <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_selected/cisd_selection.irp.f#L1>`_
+  Undocumented
+
+
+
+Needed Modules
+==============
+
+.. Do not edit this section. It was auto-generated from the
+.. NEEDED_MODULES file.
+
+* `AOs <http://github.com/LCPQ/quantum_package/tree/master/src/AOs>`_
+* `BiInts <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts>`_
+* `Bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask>`_
+* `CISD <http://github.com/LCPQ/quantum_package/tree/master/src/CISD>`_
+* `Dets <http://github.com/LCPQ/quantum_package/tree/master/src/Dets>`_
+* `Electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Electrons>`_
+* `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
+* `Hartree_Fock <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock>`_
+* `MonoInts <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts>`_
+* `MOs <http://github.com/LCPQ/quantum_package/tree/master/src/MOs>`_
+* `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_
+* `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
+* `Perturbation <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation>`_
+* `SingleRefMethod <http://github.com/LCPQ/quantum_package/tree/master/src/SingleRefMethod>`_
+* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
+* `Selectors_full <http://github.com/LCPQ/quantum_package/tree/master/src/Selectors_full>`_
+

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

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

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

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

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

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

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