Merge pull request #39 from QuantumPackage/dev-lcpq

Dev lcpq

bin/qp_set_frozen_core

@@ -13,10 +13,11 @@ zero.
 
 
 Usage:
-      qp_set_frozen_core [-q|--query] EZFIO_DIR
+      qp_set_frozen_core [-q|--query] [-l|--large] EZFIO_DIR
 
 Options:
     -q --query   Prints in the standard output the number of frozen MOs
+    -l --large   Use a large core
 
 """
 
@@ -46,7 +47,34 @@ def main(arguments):
     except:
         do_pseudo = False
 
+    large = 0
+    small = 1
+
+    size = small
+    if arguments["--large"]:
+        size = large
+
     if not do_pseudo:
+
+        if size == large:
+            for charge in ezfio.nuclei_nucl_charge:
+		if charge <= 2:
+		    pass 
+		elif charge <= 10:
+		    n_frozen += 1
+		elif charge <= 18:
+		    n_frozen += 5
+		elif charge <= 36:
+		    n_frozen += 9
+		elif charge <= 54:
+		    n_frozen += 18
+		elif charge <= 86:
+		    n_frozen += 27
+		elif charge <= 118:
+		    n_frozen += 43
+
+        if size == small:
+
             for charge in ezfio.nuclei_nucl_charge:
                 if charge < 5:
                     pass
@@ -54,8 +82,12 @@ def main(arguments):
                     n_frozen += 1
                 elif charge < 31:
                     n_frozen += 5
-            else:
+                elif charge < 49:
                     n_frozen += 9
+                elif charge < 81:
+                    n_frozen += 18
+                elif charge < 113:
+                    n_frozen += 27
 
     mo_num = ezfio.mo_basis_mo_num
 
@@ -65,10 +97,10 @@ def main(arguments):
 
     if n_frozen == 0:
         os.system("""qp_set_mo_class -a "[1-%d]" %s""" %
-                  (mo_num, sys.argv[1]))
+                  (mo_num, filename))
     else:
         os.system("""qp_set_mo_class -c "[1-%d]" -a "[%d-%d]" %s""" %
-                  (n_frozen, n_frozen+1, mo_num, sys.argv[1]))
+                  (n_frozen, n_frozen+1, mo_num, filename))
 
 
 

bin/qp_set_frozen_large_core

@@ -0,0 +1,79 @@
+#!/usr/bin/env python2
+
+
+"""
+Automatically finds n, the number of core electrons. Calls qp_set_mo_class
+setting all MOs as Active, except the n/2 first ones which are set as Core.
+If pseudo-potentials are used, all the MOs are set as Active.
+
+
+
+Usage:
+      qp_set_frozen_core [-q|--query] EZFIO_DIR
+
+Options:
+    -q --query   Prints in the standard output the number of frozen MOs
+
+"""
+
+import os
+import sys
+import os.path
+
+try:
+    import qp_path
+except ImportError:
+    print "source .quantum_package.rc"
+    raise
+
+from docopt import docopt
+from ezfio import ezfio
+
+
+def main(arguments):
+    """Main function"""
+
+    filename = arguments["EZFIO_DIR"]
+    ezfio.set_filename(filename)
+
+    n_frozen = 0
+    try:
+        do_pseudo = ezfio.pseudo_do_pseudo
+    except:
+        do_pseudo = False
+
+    if not do_pseudo:
+        for charge in ezfio.nuclei_nucl_charge:
+            if charge <= 2:
+                pass 
+            elif charge <= 10:
+                n_frozen += 1
+            elif charge <= 18:
+                n_frozen += 5
+            elif charge <= 36:
+                n_frozen += 9
+            elif charge <= 54:
+                n_frozen += 18
+            elif charge <= 86:
+                n_frozen += 27
+            elif charge <= 118:
+                n_frozen += 43
+
+    mo_num = ezfio.mo_basis_mo_num
+
+    if arguments["--query"]:
+        print n_frozen
+        sys.exit(0)
+
+    if n_frozen == 0:
+        os.system("""qp_set_mo_class -a "[1-%d]" %s""" %
+                  (mo_num, sys.argv[1]))
+    else:
+        os.system("""qp_set_mo_class -c "[1-%d]" -a "[%d-%d]" %s""" %
+                  (n_frozen, n_frozen+1, mo_num, sys.argv[1]))
+
+
+
+if __name__ == '__main__':
+    ARGUMENTS = docopt(__doc__)
+    main(ARGUMENTS)

configure

@@ -290,8 +290,7 @@ EOF
                   | sh \${QP_ROOT}/external/opam_installer.sh 
                   rm \${QP_ROOT}/external/opam_installer.sh
                   source \${OPAMROOT}/opam-init/init.sh > /dev/null 2> /dev/null || true
-                  \${QP_ROOT}/bin/opam init --disable-sandboxing --verbose \
+                  \${QP_ROOT}/bin/opam init --disable-sandboxing --verbose --yes
-                      --yes 
                   eval \$(\${QP_ROOT}/bin/opam env)
                   opam install -y \${OCAML_PACKAGES}  || exit 1
 EOF

etc/.gitignore

@@ -1 +1,2 @@
 00.qp_root
+local.rc

etc/local.rc

@@ -16,6 +16,7 @@
 # export OMP_NUM_THREADS=16
 
 # Name of the network interface to be chosen
+# export QP_NIC=lo
 # export QP_NIC=ib0
 
 

ocaml/Input_determinants_by_hand.ml

@@ -84,7 +84,7 @@ end = struct
     let n_det_old =
       Ezfio.get_determinants_n_det ()
     in
-    min n_det_old (Det_number.to_int n)
+    Det_number.to_int n
     |> Ezfio.set_determinants_n_det
   ;;
 

ocaml/MO_label.ml

@@ -6,6 +6,7 @@ type t =
 | Natural
 | Localized
 | Orthonormalized
+| MCSCF
 | None
 [@@deriving sexp]
 
@@ -16,6 +17,7 @@ let to_string = function
   | Orthonormalized -> "Orthonormalized"
   | Natural   -> "Natural"
   | Localized -> "Localized"
+  | MCSCF -> "MCSCF"
   | None      -> "None"
 ;;
 
@@ -26,7 +28,8 @@ let of_string  s =
   | "natural"   -> Natural
   | "localized" -> Localized
   | "orthonormalized" -> Orthonormalized
+  | "mcscf" -> MCSCF
   | "none"      -> None
   | _ -> (print_endline s ; failwith "MO_label should be one of:
-Guess | Orthonormalized | Canonical | Natural | Localized | None.")
+Guess | Orthonormalized | Canonical | Natural | Localized | MCSCF | None.")
 ;;

ocaml/MO_label.mli

@@ -4,6 +4,7 @@ type t =
   | Natural
   | Localized
   | Orthonormalized
+  | MCSCF
   | None
 [@@deriving sexp]
 

ocaml/element_create_db.ml

@@ -2,9 +2,6 @@ open Qptypes
 open Element
 
 let () =
-  let indices =
-    Array.init 78 (fun i -> i)
-  in
   let out_channel = 
     open_out (Qpackage.root ^ "/data/list_element.txt")
   in

scripts/generate_h_apply.py

@@ -5,6 +5,8 @@ import os
 keywords = """
 check_double_excitation
 copy_buffer
+filter_only_connected_to_hf_single
+filter_only_connected_to_hf_double
 declarations
 decls_main
 deinit_thread
@@ -205,84 +207,84 @@ class H_apply(object):
   def filter_only_2h(self):
     self["only_2h_single"] = """
 !    ! DIR$ FORCEINLINE
-     if (is_a_2h(hole).eqv. .False.) cycle
+     if (.not.is_a_2h(hole)) cycle
     """
     self["only_2h_double"] = """
 !    ! DIR$ FORCEINLINE
-     if ( is_a_2h(key).eqv. .False. )cycle
+     if (.not.is_a_2h(key))cycle
     """
 
   def filter_only_1h(self):
     self["only_1h_single"] = """
 !    ! DIR$ FORCEINLINE
-     if (is_a_1h(hole) .eqv. .False.) cycle
+     if (.not.is_a_1h(hole)) cycle
     """
     self["only_1h_double"] = """
 !    ! DIR$ FORCEINLINE
-     if (is_a_1h(key) .eqv. .False.) cycle
+     if (.not.is_a_1h(key)  ) cycle
     """
 
   def filter_only_1p(self):
     self["only_1p_single"] = """
 !    ! DIR$ FORCEINLINE
-     if ( is_a_1p(hole) .eqv. .False.) cycle
+     if (.not. is_a_1p(hole) ) cycle
     """
     self["only_1p_double"] = """
 !    ! DIR$ FORCEINLINE
-     if ( is_a_1p(key) .eqv. .False.) cycle
+     if (.not. is_a_1p(key) ) cycle
     """
 
   def filter_only_2h1p(self):
     self["only_2h1p_single"] = """
 !    ! DIR$ FORCEINLINE
-     if ( is_a_2h1p(hole) .eqv. .False.) cycle
+     if (.not. is_a_2h1p(hole) ) cycle
     """
     self["only_2h1p_double"] = """
 !    ! DIR$ FORCEINLINE
-     if (is_a_2h1p(key) .eqv. .False.) cycle
+     if (.not.is_a_2h1p(key) ) cycle
     """
 
 
   def filter_only_2p(self):
     self["only_2p_single"] = """
 !    ! DIR$ FORCEINLINE
-     if (is_a_2p(hole).eqv. .False.) cycle
+     if (.not.is_a_2p(hole)) cycle
     """
     self["only_2p_double"] = """
 !    ! DIR$ FORCEINLINE
-     if (is_a_2p(key).eqv. .False.) cycle
+     if (.not.is_a_2p(key)) cycle
     """
 
 
   def filter_only_1h1p(self):
     self["filter_only_1h1p_single"] = """
-!    ! DIR$ FORCEINLINE
+     if (.not.is_a_1h1p(hole)) cycle
-     if (is_a_1h1p(hole).eqv..False.) cycle
     """
     self["filter_only_1h1p_double"] = """
-!    ! DIR$ FORCEINLINE
+     if (.not.is_a_1h1p(key)) cycle
-     if (is_a_1h1p(key).eqv..False.) cycle
     """
 
+
+
   def filter_only_2h2p(self):
     self["filter_only_2h2p_single"] = """
 !    ! DIR$ FORCEINLINE
-     if (is_a_two_holes_two_particles(hole).eqv..False.) cycle
+     if (.not.is_a_two_holes_two_particles(hole)) cycle
     """
     self["filter_only_2h2p_double"] = """
 !    ! DIR$ FORCEINLINE
-     if (is_a_two_holes_two_particles(key).eqv..False.) cycle
+     if (.not.is_a_two_holes_two_particles(key)) cycle
     """
 
 
   def filter_only_1h2p(self):
     self["filter_only_1h2p_single"] = """
 !    ! DIR$ FORCEINLINE
-     if (is_a_1h2p(hole).eqv..False.) cycle
+     if (.not.is_a_1h2p(hole)) cycle
     """
     self["filter_only_1h2p_double"] = """
 !    ! DIR$ FORCEINLINE
-     if (is_a_1h2p(key).eqv..False.) cycle
+     if (.not.is_a_1h2p(key)) cycle
     """
 
 
@@ -294,6 +296,16 @@ class H_apply(object):
      if (is_a_two_holes_two_particles(hole)) cycle
     """
 
+  def filter_only_connected_to_hf(self):
+    self["filter_only_connected_to_hf_single"] = """
+     call connected_to_hf(hole,yes_no)
+     if (.not.yes_no) cycle
+    """
+    self["filter_only_connected_to_hf_double"] = """
+     call connected_to_hf(key,yes_no)
+     if (.not.yes_no) cycle
+    """
+
 
   def set_perturbation(self,pert):
     if self.perturbation is not None:

src/ao_two_e_ints/map_integrals.irp.f

@@ -19,6 +19,10 @@ END_PROVIDER
 subroutine two_e_integrals_index(i,j,k,l,i1)
   use map_module
   implicit none
+  BEGIN_DOC
+! Gives a unique index for i,j,k,l using permtuation symmetry.
+! i <-> k, j <-> l, and (i,k) <-> (j,l)
+  END_DOC
   integer, intent(in)            :: i,j,k,l
   integer(key_kind), intent(out) :: i1
   integer(key_kind)              :: p,q,r,s,i2
@@ -36,14 +40,25 @@ end
 subroutine two_e_integrals_index_reverse(i,j,k,l,i1)
   use map_module
   implicit none
+  BEGIN_DOC
+! Computes the 4 indices $i,j,k,l$ from a unique index $i_1$.
+! For 2 indices $i,j$ and $i \le j$, we have
+! $p = i(i-1)/2 + j$.
+! The key point is that because $j < i$,
+! $i(i-1)/2 < p \le i(i+1)/2$. So $i$ can be found by solving
+! $i^2 - i - 2p=0$. One obtains $i=1 + \sqrt{1+8p}/2$
+! and $j = p - i(i-1)/2$.
+! This rule is applied 3 times. First for the symmetry of the
+! pairs (i,k) and (j,l), and then for the symmetry within each pair.
+  END_DOC
   integer, intent(out)           :: i(8),j(8),k(8),l(8)
   integer(key_kind), intent(in)  :: i1
   integer(key_kind)              :: i2,i3
   i = 0
-  i2   = ceiling(0.5d0*(dsqrt(8.d0*dble(i1)+1.d0)-1.d0))
+  i2   = ceiling(0.5d0*(dsqrt(dble(shiftl(i1,3)+1))-1.d0))
-  l(1) = ceiling(0.5d0*(dsqrt(8.d0*dble(i2)+1.d0)-1.d0))
+  l(1) = ceiling(0.5d0*(dsqrt(dble(shiftl(i2,3)+1))-1.d0))
   i3   = i1 - shiftr(i2*i2-i2,1)
-  k(1) = ceiling(0.5d0*(dsqrt(8.d0*dble(i3)+1.d0)-1.d0))
+  k(1) = ceiling(0.5d0*(dsqrt(dble(shiftl(i3,3)+1))-1.d0))
   j(1) = int(i2 - shiftr(l(1)*l(1)-l(1),1),4)
   i(1) = int(i3 - shiftr(k(1)*k(1)-k(1),1),4)
 
@@ -95,16 +110,18 @@ subroutine two_e_integrals_index_reverse(i,j,k,l,i1)
       endif
     enddo
   enddo
-  do ii=1,8
+! This has been tested with up to 1000 AOs, and all the reverse indices are
-    if (i(ii) /= 0) then
+! correct ! We can remove the test
-      call two_e_integrals_index(i(ii),j(ii),k(ii),l(ii),i2)
+!    do ii=1,8
-      if (i1 /= i2) then
+!      if (i(ii) /= 0) then
-        print *,  i1, i2
+!        call two_e_integrals_index(i(ii),j(ii),k(ii),l(ii),i2)
-        print *,  i(ii), j(ii), k(ii), l(ii)
+!        if (i1 /= i2) then
-        stop 'two_e_integrals_index_reverse failed'
+!          print *,  i1, i2
-      endif
+!          print *,  i(ii), j(ii), k(ii), l(ii)
-    endif
+!          stop 'two_e_integrals_index_reverse failed'
-  enddo
+!        endif
+!      endif
+!    enddo
 
 
 end
@@ -196,6 +213,7 @@ subroutine get_ao_two_e_integrals(j,k,l,sze,out_val)
   BEGIN_DOC
   ! Gets multiple AO bi-electronic integral from the AO map .
   ! All i are retrieved for j,k,l fixed.
+  ! physicist convention : <ij|kl> 
   END_DOC
   implicit none
   integer, intent(in)            :: j,k,l, sze

src/aux_quantities/EZFIO.cfg

@@ -24,3 +24,17 @@ type: double precision
 size: (mo_basis.mo_num,mo_basis.mo_num,determinants.n_states)
 
 
+[data_one_e_dm_alpha_ao]
+interface: ezfio, provider
+doc: Alpha one body density matrix on the |AO| basis computed with the wave function
+type: double precision 
+size: (ao_basis.ao_num,ao_basis.ao_num,determinants.n_states)
+
+
+[data_one_e_dm_beta_ao]
+interface: ezfio, provider
+doc: Beta one body density matrix on the |AO| basis computed with the wave function
+type: double precision 
+size: (ao_basis.ao_num,ao_basis.ao_num,determinants.n_states)
+
+

src/bitmask/bitmasks.irp.f

@@ -466,35 +466,17 @@ END_PROVIDER
 
 
  BEGIN_PROVIDER [integer(bit_kind), reunion_of_core_inact_act_bitmask, (N_int,2)]
-&BEGIN_PROVIDER [ integer, n_core_inact_act_orb ]
  implicit none
  BEGIN_DOC
  ! Reunion of the core, inactive and active bitmasks
  END_DOC
  integer :: i,j
 
- n_core_inact_act_orb = 0
  do i = 1, N_int
-  reunion_of_core_inact_act_bitmask(i,1) = ior(reunion_of_core_inact_bitmask(i,1),cas_bitmask(i,1,1))
+  reunion_of_core_inact_act_bitmask(i,1) = ior(reunion_of_core_inact_bitmask(i,1),act_bitmask(i,1))
-  reunion_of_core_inact_act_bitmask(i,2) = ior(reunion_of_core_inact_bitmask(i,2),cas_bitmask(i,2,1))
+  reunion_of_core_inact_act_bitmask(i,2) = ior(reunion_of_core_inact_bitmask(i,2),act_bitmask(i,2))
-  n_core_inact_act_orb +=popcnt(reunion_of_core_inact_act_bitmask(i,1))
  enddo
  END_PROVIDER
- BEGIN_PROVIDER [ integer, list_core_inact_act, (n_core_inact_act_orb)]
-&BEGIN_PROVIDER [ integer, list_core_inact_act_reverse, (mo_num)]
-  implicit none
- integer :: occ_inact(N_int*bit_kind_size)
- integer :: itest,i
- occ_inact = 0
- call bitstring_to_list(reunion_of_core_inact_act_bitmask(1,1), occ_inact(1), itest, N_int)
- list_inact_reverse = 0
- do i = 1, n_core_inact_act_orb
-  list_core_inact_act(i) = occ_inact(i)
-  list_core_inact_act_reverse(occ_inact(i)) = i
- enddo
-END_PROVIDER
-
-
 
 
  BEGIN_PROVIDER [ integer(bit_kind), reunion_of_bitmask, (N_int,2)]
@@ -563,8 +545,8 @@ END_PROVIDER
  END_DOC
  integer :: i,j
  do i = 1, N_int
-  reunion_of_cas_inact_bitmask(i,1) = ior(cas_bitmask(i,1,1),inact_bitmask(i,1))
+  reunion_of_cas_inact_bitmask(i,1) = ior(act_bitmask(i,1),inact_bitmask(i,1))
-  reunion_of_cas_inact_bitmask(i,2) = ior(cas_bitmask(i,2,1),inact_bitmask(i,2))
+  reunion_of_cas_inact_bitmask(i,2) = ior(act_bitmask(i,2),inact_bitmask(i,2))
  enddo
  END_PROVIDER
 

src/bitmask/core_inact_act_virt.irp.f

@@ -194,3 +194,53 @@ END_PROVIDER
 
 END_PROVIDER 
 
+BEGIN_PROVIDER [integer, n_inact_act_orb ]
+ implicit none
+ n_inact_act_orb = (n_inact_orb+n_act_orb)
+
+END_PROVIDER 
+
+BEGIN_PROVIDER [integer, list_inact_act, (n_inact_act_orb)]
+ integer :: i,itmp
+ itmp = 0
+ do i = 1, n_inact_orb
+  itmp += 1
+  list_inact_act(itmp) = list_inact(i) 
+ enddo
+ do i = 1, n_act_orb
+  itmp += 1
+  list_inact_act(itmp) = list_act(i) 
+ enddo
+END_PROVIDER 
+
+BEGIN_PROVIDER [integer, n_core_inact_act_orb ]
+ implicit none
+ n_core_inact_act_orb = (n_core_orb + n_inact_orb + n_act_orb)
+
+END_PROVIDER 
+
+ BEGIN_PROVIDER [integer, list_core_inact_act, (n_core_inact_act_orb)]
+&BEGIN_PROVIDER [ integer, list_core_inact_act_reverse, (n_core_inact_act_orb)]
+ integer :: i,itmp
+ itmp = 0
+ do i = 1, n_core_orb
+  itmp += 1
+  list_core_inact_act(itmp) = list_core(i) 
+ enddo
+ do i = 1, n_inact_orb
+  itmp += 1
+  list_core_inact_act(itmp) = list_inact(i) 
+ enddo
+ do i = 1, n_act_orb
+  itmp += 1
+  list_core_inact_act(itmp) = list_act(i) 
+ enddo
+
+ integer :: occ_inact(N_int*bit_kind_size)
+ occ_inact = 0
+ call bitstring_to_list(reunion_of_core_inact_act_bitmask(1,1), occ_inact(1), itest, N_int)
+ list_inact_reverse = 0
+ do i = 1, n_core_inact_act_orb
+  list_core_inact_act_reverse(occ_inact(i)) = i
+ enddo
+END_PROVIDER 

src/cipsi/pt2_stoch_routines.irp.f

@@ -333,6 +333,14 @@ subroutine ZMQ_pt2(E, pt2,relative_error, error, variance, norm, N_in)
     pt2(k) = 0.d0
   enddo
 
+  ! Adjust PT2 weights for next selection
+  double precision :: pt2_avg
+  pt2_avg = sum(pt2) / dble(N_states)
+  do k=1,N_states
+    pt2_match_weight(k) *= (pt2(k)/pt2_avg)**2
+  enddo
+  SOFT_TOUCH pt2_match_weight
+
 end subroutine
 
 

src/cipsi/selection.irp.f

@@ -1,11 +1,20 @@
 use bitmasks
 
+BEGIN_PROVIDER [ double precision, pt2_match_weight, (N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! Weights adjusted along the selection to make the PT2 contributions
+ ! of each state coincide.
+ END_DOC
+ pt2_match_weight = 1.d0
+END_PROVIDER
+
 BEGIN_PROVIDER [ double precision, selection_weight, (N_states) ]
    implicit none
    BEGIN_DOC
    ! Weights used in the selection criterion
    END_DOC
-   selection_weight(1:N_states) = c0_weight(1:N_states)
+   selection_weight(1:N_states) = c0_weight(1:N_states) * pt2_match_weight(1:N_states)
 END_PROVIDER
 
 
@@ -618,6 +627,11 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
           sum_e_pert = sum_e_pert + e_pert * selection_weight(istate)
 !        endif
       end do
+      if(pseudo_sym)then
+       if(dabs(mat(1, p1, p2)).lt.thresh_sym)then 
+        sum_e_pert = 10.d0
+       endif
+      endif
 
       if(sum_e_pert <= buf%mini) then
         call add_to_selection_buffer(buf, det, sum_e_pert)

src/cipsi/zmq_selection.irp.f

@@ -131,6 +131,14 @@ subroutine ZMQ_selection(N_in, pt2, variance, norm)
     norm(k) = norm(k) * f(k)
   enddo
 
+  ! Adjust PT2 weights for next selection
+  double precision :: pt2_avg
+  pt2_avg = sum(pt2) / dble(N_states)
+  do k=1,N_states
+    pt2_match_weight(k) *= (pt2(k)/pt2_avg)**2
+  enddo
+  SOFT_TOUCH pt2_match_weight
+
 end subroutine
 
 

src/cis/cis.irp.f

@@ -57,7 +57,11 @@ subroutine run
   implicit none
   integer                        :: i
 
+  if(pseudo_sym)then
+   call H_apply_cis_sym
+  else
    call H_apply_cis
+  endif
   print *,  'N_det = ', N_det
   print*,'******************************'
   print *,  'Energies  of the states:'

src/cis/h_apply.irp.f

@@ -5,5 +5,10 @@ BEGIN_SHELL [ /usr/bin/env python2 ]
 from generate_h_apply import H_apply
 H = H_apply("cis",do_double_exc=False)
 print H
+
+H = H_apply("cis_sym",do_double_exc=False)
+H.filter_only_connected_to_hf()
+print H
+
 END_SHELL
 

src/cisd/cisd.irp.f

@@ -53,7 +53,11 @@ subroutine run
   implicit none
   integer :: i
 
+  if(pseudo_sym)then
+   call H_apply_cisd_sym
+  else
    call H_apply_cisd
+  endif
   print *,  'N_det = ', N_det
   print*,'******************************'
   print *,  'Energies  of the states:'

src/cisd/h_apply.irp.f

@@ -5,5 +5,9 @@ BEGIN_SHELL [ /usr/bin/env python2 ]
 from generate_h_apply import H_apply
 H = H_apply("cisd",do_double_exc=True)
 print H
+
+H = H_apply("cisd_sym",do_double_exc=True)
+H.filter_only_connected_to_hf()
+print H
 END_SHELL
 

src/davidson/diagonalization_hs2_dressed.irp.f

@@ -157,7 +157,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: U, W, S, y, y_s, S_d, h, lambda
   if (N_st_diag*3 > sze) then
     print *,  'error in Davidson :'
-    print *,  'Increase n_det_max_jacobi to ', N_st_diag*3
+    print *,  'Increase n_det_max_full to ', N_st_diag*3
     stop -1
   endif
 

src/density_for_dft/density_for_dft.irp.f

@@ -9,6 +9,8 @@ BEGIN_PROVIDER [double precision, one_e_dm_mo_alpha_for_dft, (mo_num,mo_num, N_s
   one_e_dm_mo_alpha_for_dft = data_one_e_dm_alpha_mo + damping_for_rs_dft * delta_alpha
  else if (density_for_dft .EQ. "input_density")then
   one_e_dm_mo_alpha_for_dft = data_one_e_dm_alpha_mo
+ else if (density_for_dft .EQ. "input_density_ao")then
+  call ao_to_mo(data_one_e_dm_alpha_mo,size(data_one_e_dm_alpha_mo,1),one_e_dm_mo_alpha_for_dft,size(one_e_dm_mo_alpha_for_dft,1))
  else if (density_for_dft .EQ. "WFT")then
   provide mo_coef
   one_e_dm_mo_alpha_for_dft = one_e_dm_mo_alpha
@@ -58,6 +60,8 @@ BEGIN_PROVIDER [double precision, one_e_dm_mo_beta_for_dft, (mo_num,mo_num, N_st
   one_e_dm_mo_beta_for_dft = data_one_e_dm_beta_mo + damping_for_rs_dft * delta_beta
  else if (density_for_dft .EQ. "input_density")then
   one_e_dm_mo_beta_for_dft = data_one_e_dm_beta_mo
+ else if (density_for_dft .EQ. "input_density_ao")then
+  call ao_to_mo(data_one_e_dm_beta_mo,size(data_one_e_dm_beta_mo,1),one_e_dm_mo_beta_for_dft,size(one_e_dm_mo_beta_for_dft,1))
  else if (density_for_dft .EQ. "WFT")then
   provide mo_coef
   one_e_dm_mo_beta_for_dft = one_e_dm_mo_beta
@@ -119,6 +123,11 @@ END_PROVIDER
 
  one_e_dm_alpha_ao_for_dft = 0.d0
  one_e_dm_beta_ao_for_dft = 0.d0
+
+ if (density_for_dft .EQ. "input_density_ao")then
+  one_e_dm_alpha_ao_for_dft = data_one_e_dm_alpha_ao 
+  one_e_dm_beta_ao_for_dft = data_one_e_dm_beta_ao 
+ else
   do istate = 1, N_states
    call mo_to_ao_no_overlap( one_e_dm_mo_alpha_for_dft(1,1,istate), &
                              size(one_e_dm_mo_alpha_for_dft,1),     &
@@ -129,6 +138,7 @@ END_PROVIDER
                              one_e_dm_beta_ao_for_dft(1,1,istate), &
                              size(one_e_dm_beta_ao_for_dft,1) )
   enddo
+ endif
 
 END_PROVIDER
 

src/determinants/EZFIO.cfg

@@ -89,3 +89,16 @@ doc: Weight of the states in state-average calculations.
 interface: ezfio
 size: (determinants.n_states)
 
+
+[thresh_sym]
+type: Threshold
+doc: Thresholds to check if a determinant is connected with HF 
+interface: ezfio,provider,ocaml
+default: 1.e-15
+
+[pseudo_sym]
+type: logical
+doc: If |true|, discard any Slater determinants with an interaction smaller than thresh_sym with HF. 
+interface: ezfio,provider,ocaml
+default: False
+

src/determinants/h_apply.template.f

@@ -150,6 +150,7 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl
   logical :: is_a_2h1p
   logical :: is_a_2h
   logical :: b_cycle
+  logical :: yes_no
   check_double_excitation = .True.
   iproc = iproc_in
 
@@ -284,6 +285,7 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl
           $only_1h_double
           $only_1p_double
           $only_2h1p_double
+          $filter_only_connected_to_hf_double 
           key_idx += 1
           do k=1,N_int
             keys_out(k,1,key_idx) = key(k,1)
@@ -339,6 +341,7 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl
         $only_1h_double
         $only_1p_double
         $only_2h1p_double
+        $filter_only_connected_to_hf_double
         key_idx += 1
         do k=1,N_int
           keys_out(k,1,key_idx) = key(k,1)
@@ -412,6 +415,7 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,fock_diag_tmp,i_generato
   logical :: is_a_1h
   logical :: is_a_1p
   logical :: is_a_2p
+  logical :: yes_no
 
   do k=1,N_int
     key_mask(k,1) = 0_bit_kind
@@ -493,6 +497,7 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,fock_diag_tmp,i_generato
       $filter_only_1h1p_single
       $filter_only_1h2p_single
       $filter_only_2h2p_single
+      $filter_only_connected_to_hf_single 
       key_idx += 1
       do k=1,N_int
         keys_out(k,1,key_idx) = hole(k,1)

src/determinants/slater_rules.irp.f

@@ -2257,3 +2257,38 @@ subroutine i_H_j_double_alpha_beta(key_i,key_j,Nint,hij)
 end
 
 
+subroutine connected_to_hf(key_i,yes_no)
+ implicit none 
+ use bitmasks
+ integer(bit_kind), intent(in)  :: key_i(N_int,2)
+ logical , intent(out) :: yes_no
+ double precision :: hij,phase
+ integer          :: exc(0:2,2,2)
+ integer          :: degree
+ integer          :: m,p
+ yes_no = .True.
+ call get_excitation_degree(ref_bitmask,key_i,degree,N_int)
+ if(degree == 2)then
+  call i_H_j(ref_bitmask,key_i,N_int,hij)
+  if(dabs(hij) .lt. thresh_sym)then
+   yes_no = .False. 
+  endif
+ else if(degree == 1)then  
+  call get_single_excitation(ref_bitmask,key_i,exc,phase,N_int)
+  ! Single alpha
+  if (exc(0,1,1) == 1) then
+    m = exc(1,1,1)
+    p = exc(1,2,1)
+  ! Single beta
+  else
+    m = exc(1,1,2)
+    p = exc(1,2,2)
+  endif
+  hij = mo_one_e_integrals(m,p)
+  if(dabs(hij) .lt. thresh_sym)then
+   yes_no = .False. 
+  endif
+ else if(degree == 0)then
+  yes_no = .True.
+ endif
+end

src/tools/rotate_mos.irp.f

@@ -0,0 +1,18 @@
+program rotate_mos
+ implicit none
+ integer :: iorb,jorb
+ read(5,*)iorb,jorb
+ double precision, allocatable :: mo_coef_tmp(:,:)
+ allocate(mo_coef_tmp(ao_num,mo_num))
+ mo_coef_tmp = mo_coef
+ integer :: i,j
+ double precision :: dsqrt2_inv
+ dsqrt2_inv = 1.d0/dsqrt(2.d0)
+ do i = 1, ao_num
+  mo_coef(i,iorb) = dsqrt2_inv * ( mo_coef_tmp(i,iorb) + mo_coef_tmp(i,jorb) )
+  mo_coef(i,jorb) = dsqrt2_inv * ( mo_coef_tmp(i,iorb) - mo_coef_tmp(i,jorb) )
+ enddo
+ touch mo_coef 
+ call save_mos
+
+end

src/tools/save_one_e_dm.irp.f

@@ -1,15 +1,15 @@
 program save_one_e_dm
   implicit none
  BEGIN_DOC
-! Program that computes the one body density on the |MO| basis
+! Program that computes the one body density on the |MO| and |AO| basis
 ! for $\alpha$ and $\beta$ electrons from the wave function
 ! stored in the |EZFIO| directory, and then saves it into the
 ! :ref:`module_aux_quantities`.
 !
 ! Then, the global variable :option:`aux_quantities data_one_e_dm_alpha_mo`
-! and :option:`aux_quantities data_one_e_dm_beta_mo` will automatically
+! and :option:`aux_quantities data_one_e_dm_beta_mo` (and the corresponding for |AO|)
-! read this density in the next calculation. This can be used to perform
+! will automatically ! read this density in the next calculation. 
-! damping on the density in |RSDFT| calculations (see
+! This can be used to perform damping on the density in |RSDFT| calculations (see
 ! :ref:`module_density_for_dft`).
  END_DOC
   read_wf = .True.
@@ -25,4 +25,6 @@ subroutine routine_save_one_e_dm
  END_DOC
  call ezfio_set_aux_quantities_data_one_e_dm_alpha_mo(one_e_dm_mo_alpha)
  call ezfio_set_aux_quantities_data_one_e_dm_beta_mo(one_e_dm_mo_beta)
+ call ezfio_set_aux_quantities_data_one_e_dm_alpha_ao(one_e_dm_ao_alpha)
+ call ezfio_set_aux_quantities_data_one_e_dm_beta_ao(one_e_dm_ao_beta)
 end