Merge pull request #26 from LCPQ/master

Merging

ocaml/qptypes_generator.ml

@@ -150,12 +150,12 @@ end = struct
 
   let to_string = function
   | Huckel -> \"Huckel\"
-  | HCore  -> \"Hcore\"
+  | HCore  -> \"HCore\"
 
   let of_string  s = 
-    match s with
+    match (String.lowercase s) with
-    | \"Huckel\" -> Huckel
+    | \"huckel\" -> Huckel
-    | \"Hcore\"  -> HCore
+    | \"hcore\"  -> HCore
     | _ -> failwith (\"Wrong Guess type : \"^s)
 
 end
@@ -176,10 +176,10 @@ end = struct
   | Write  -> \"Write\"
   | None   -> \"None\"
   let of_string  s = 
-    match s with
+    match (String.lowercase s) with
-    | \"Read\"  -> Read
+    | \"read\"  -> Read
-    | \"Write\" -> Write
+    | \"write\" -> Write
-    | \"None\"  -> None
+    | \"none\"  -> None
     | _ -> failwith (\"Wrong IO type : \"^s)
 
 end

plugins/Full_CI/.gitignore

@@ -22,6 +22,7 @@ Properties
 Pseudo
 Selectors_full
 Utils
+exc_degree
 ezfio_interface.irp.f
 full_ci
 full_ci_no_skip

plugins/MRCC/H_apply.irp.f

@@ -1,28 +0,0 @@
-use bitmasks
-BEGIN_SHELL [ /usr/bin/env python ]
-from generate_h_apply import *
-
-s = H_apply("mrcc")
-s.data["parameters"] = ", delta_ij_, delta_ii_,Ndet_cas, Ndet_non_cas"
-s.data["declarations"] += """
-    integer, intent(in) :: Ndet_cas,Ndet_non_cas
-    double precision, intent(in) :: delta_ij_(Ndet_cas,Ndet_non_cas,*)
-    double precision, intent(in) :: delta_ii_(Ndet_cas,*)
-"""
-s.data["keys_work"] = "call mrcc_dress(delta_ij_,delta_ii_,Ndet_cas,Ndet_non_cas,i_generator,key_idx,keys_out,N_int,iproc)"
-s.data["params_post"] += ", delta_ij_, delta_ii_, Ndet_cas, Ndet_non_cas"
-s.data["params_main"] += "delta_ij_, delta_ii_, Ndet_cas, Ndet_non_cas"
-s.data["decls_main"] += """
-    integer, intent(in) :: Ndet_cas,Ndet_non_cas
-    double precision, intent(in) :: delta_ij_(Ndet_cas,Ndet_non_cas,*)
-    double precision, intent(in) :: delta_ii_(Ndet_cas,*)
-"""
-s.data["finalization"] = ""
-s.data["copy_buffer"] = ""
-s.data["generate_psi_guess"] = ""
-s.data["size_max"] = "3072"
-print s
-
-
-END_SHELL
-

plugins/MRCC/NEEDED_CHILDREN_MODULES

@@ -1 +0,0 @@
-Perturbation Selectors_full Generators_full

plugins/MRCC/mrcc.irp.f

@@ -1,70 +0,0 @@
-program mrcc
-  implicit none
-  read_wf = .True.
-  TOUCH read_wf
-  call run
-  call run_mrcc
-!  call run_mrcc_test
-end
-
-subroutine run
-  implicit none
-
-  integer :: i,j
-  print *,  'CAS'
-  print *,  '==='
-  do i=1,N_det_cas
-    print *,  psi_cas_coef(i,:)
-    call debug_det(psi_cas(1,1,i),N_int)
-  enddo
-
-!  print *,  'SD'
-!  print *,  '=='
-!  do i=1,N_det_non_cas
-!    print *,  psi_non_cas_coef(i,:)
-!    call debug_det(psi_non_cas(1,1,i),N_int)
-!  enddo
-  call write_double(6,ci_energy(1),"Initial CI energy")
-end
-subroutine run_mrcc_test
-  implicit none
-  integer :: i,j
-  double precision :: pt2
-  pt2 = 0.d0
-  do j=1,N_det
-   do i=1,N_det
-     pt2 += psi_coef(i,1)*psi_coef(j,1) * delta_ij(i,j,1)
-   enddo
-  enddo
-  print *,  ci_energy(1)
-  print *,  ci_energy(1)+pt2
-end
-subroutine run_mrcc
-  implicit none
-  integer :: i,j
-
-  double precision :: E_new, E_old, delta_e
-  integer :: iteration
-  E_new = 0.d0
-  delta_E = 1.d0
-  iteration = 0
-  do while (delta_E > 1.d-10)
-    iteration += 1
-    print *,  '===========================' 
-    print *,  'MRCC Iteration', iteration
-    print *,  '===========================' 
-    print *,  ''
-    E_old = sum(ci_energy_dressed)
-    call write_double(6,ci_energy_dressed(1),"MRCC energy")
-    call diagonalize_ci_dressed
-    E_new = sum(ci_energy_dressed)
-    delta_E = dabs(E_new - E_old)
-    if (iteration > 20) then
-      exit
-    endif
-  enddo
-  call write_double(6,ci_energy_dressed(1),"Final MRCC energy")
-  call ezfio_set_mrcc_energy(ci_energy_dressed(1))
-!  call save_wavefunction
-
-end

plugins/MRCC_CASSD/.gitignore

@@ -0,0 +1,32 @@
+# Automatically created by $QP_ROOT/scripts/module/module_handler.py 
+.ninja_deps
+.ninja_log
+AO_Basis
+Bitmask
+Determinants
+Electrons
+Ezfio_files
+Generators_full
+Hartree_Fock
+IRPF90_man
+IRPF90_temp
+Integrals_Bielec
+Integrals_Monoelec
+MOGuess
+MO_Basis
+MRCC_Utils
+Makefile
+Makefile.depend
+Nuclei
+Perturbation
+Properties
+Pseudo
+Psiref_CAS
+Psiref_Utils
+Selectors_full
+Utils
+ezfio_interface.irp.f
+irpf90.make
+irpf90_entities
+mrcc_cassd
+tags

plugins/MRCC_CASSD/EZFIO.cfg

@@ -0,0 +1,4 @@
+[energy]
+type: double precision
+doc: Calculated energy
+interface: ezfio

plugins/MRCC_CASSD/NEEDED_CHILDREN_MODULES

@@ -0,0 +1 @@
+Perturbation Selectors_full Generators_full Psiref_CAS MRCC_Utils

plugins/MRCC_CASSD/README.rst

@@ -0,0 +1,33 @@
+===========
+MRCC Module
+===========
+
+MRCC as a coupled cluster on a CAS+SD wave function.
+
+Needed Modules
+==============
+
+.. Do not edit this section. It was auto-generated from the
+.. by the `update_README.py` script.
+
+.. image:: tree_dependency.png
+
+* `Perturbation <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation>`_
+* `Selectors_full <http://github.com/LCPQ/quantum_package/tree/master/src/Selectors_full>`_
+* `Generators_full <http://github.com/LCPQ/quantum_package/tree/master/src/Generators_full>`_
+* `Psiref_CAS <http://github.com/LCPQ/quantum_package/tree/master/src/Psiref_CAS>`_
+* `MRCC_Utils <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils>`_
+
+Documentation
+=============
+
+.. Do not edit this section. It was auto-generated from the
+.. by the `update_README.py` script.
+
+`mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_CASSD/mrcc_cassd.irp.f#L1>`_
+  Undocumented
+
+
+`print_cas_coefs <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_CASSD/mrcc_cassd.irp.f#L11>`_
+  Undocumented
+

plugins/MRCC_CASSD/mrcc_cassd.irp.f

@@ -0,0 +1,24 @@
+program mrcc
+  implicit none
+  if (.not.read_wf) then
+    print *,  'read_wf has to be true.'
+    stop 1
+  endif
+  call print_cas_coefs
+  call run_mrcc
+end
+
+subroutine print_cas_coefs
+  implicit none
+
+  integer :: i,j
+  print *,  'CAS'
+  print *,  '==='
+  do i=1,N_det_cas
+    print *,  psi_cas_coef(i,:)
+    call debug_det(psi_cas(1,1,i),N_int)
+  enddo
+
+  call write_double(6,ci_energy(1),"Initial CI energy")
+end
+

plugins/MRCC_Utils/.gitignore

@@ -0,0 +1,31 @@
+# Automatically created by $QP_ROOT/scripts/module/module_handler.py 
+.ninja_deps
+.ninja_log
+AO_Basis
+Bitmask
+Determinants
+Electrons
+Ezfio_files
+Generators_full
+Hartree_Fock
+IRPF90_man
+IRPF90_temp
+Integrals_Bielec
+Integrals_Monoelec
+MOGuess
+MO_Basis
+Makefile
+Makefile.depend
+Nuclei
+Perturbation
+Properties
+Pseudo
+Psiref_CAS
+Psiref_Utils
+Selectors_full
+Utils
+ezfio_interface.irp.f
+irpf90.make
+irpf90_entities
+mrcc_general
+tags

plugins/MRCC_Utils/H_apply.irp.f

@@ -0,0 +1,28 @@
+use bitmasks
+BEGIN_SHELL [ /usr/bin/env python ]
+from generate_h_apply import *
+
+s = H_apply("mrcc")
+s.data["parameters"] = ", delta_ij_, delta_ii_,Ndet_ref, Ndet_non_ref"
+s.data["declarations"] += """
+    integer, intent(in) :: Ndet_ref,Ndet_non_ref
+    double precision, intent(in) :: delta_ij_(Ndet_ref,Ndet_non_ref,*)
+    double precision, intent(in) :: delta_ii_(Ndet_ref,*)
+"""
+s.data["keys_work"] = "call mrcc_dress(delta_ij_,delta_ii_,Ndet_ref,Ndet_non_ref,i_generator,key_idx,keys_out,N_int,iproc)"
+s.data["params_post"] += ", delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref"
+s.data["params_main"] += "delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref"
+s.data["decls_main"] += """
+    integer, intent(in) :: Ndet_ref,Ndet_non_ref
+    double precision, intent(in) :: delta_ij_(Ndet_ref,Ndet_non_ref,*)
+    double precision, intent(in) :: delta_ii_(Ndet_ref,*)
+"""
+s.data["finalization"] = ""
+s.data["copy_buffer"] = ""
+s.data["generate_psi_guess"] = ""
+s.data["size_max"] = "3072"
+print s
+
+
+END_SHELL
+

plugins/MRCC_Utils/NEEDED_CHILDREN_MODULES

@@ -0,0 +1 @@
+Perturbation Selectors_full Generators_full Psiref_Utils

plugins/MRCC_Utils/README.rst

@@ -13,6 +13,7 @@ Needed Modules
 * `Perturbation <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation>`_
 * `Selectors_full <http://github.com/LCPQ/quantum_package/tree/master/src/Selectors_full>`_
 * `Generators_full <http://github.com/LCPQ/quantum_package/tree/master/src/Generators_full>`_
+* `Psiref_Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Psiref_Utils>`_
 
 Documentation
 =============
@@ -20,23 +21,23 @@ Documentation
 .. Do not edit this section. It was auto-generated from the
 .. by the `update_README.py` script.
 
-`ci_eigenvectors_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_utils.irp.f#L76>`_
+`ci_eigenvectors_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L76>`_
   Eigenvectors/values of the CI matrix
 
 
-`ci_eigenvectors_s2_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_utils.irp.f#L77>`_
+`ci_eigenvectors_s2_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L77>`_
   Eigenvectors/values of the CI matrix
 
 
-`ci_electronic_energy_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_utils.irp.f#L75>`_
+`ci_electronic_energy_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L75>`_
   Eigenvectors/values of the CI matrix
 
 
-`ci_energy_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_utils.irp.f#L142>`_
+`ci_energy_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L142>`_
   N_states lowest eigenvalues of the dressed CI matrix
 
 
-`davidson_diag_hjj_mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/davidson.irp.f#L56>`_
+`davidson_diag_hjj_mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/davidson.irp.f#L56>`_
   Davidson diagonalization with specific diagonal elements of the H matrix
   .br
   H_jj : specific diagonal H matrix elements to diagonalize de Davidson
@@ -57,7 +58,7 @@ Documentation
   Initial guess vectors are not necessarily orthonormal
 
 
-`davidson_diag_mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/davidson.irp.f#L1>`_
+`davidson_diag_mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/davidson.irp.f#L1>`_
   Davidson diagonalization.
   .br
   dets_in : bitmasks corresponding to determinants
@@ -76,45 +77,45 @@ Documentation
   Initial guess vectors are not necessarily orthonormal
 
 
-`delta_ii <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_utils.irp.f#L39>`_
+`delta_ii <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L39>`_
   Dressing matrix in N_det basis
 
 
-`delta_ij <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_utils.irp.f#L38>`_
+`delta_ij <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L38>`_
   Dressing matrix in N_det basis
 
 
-`diagonalize_ci_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_utils.irp.f#L157>`_
+`diagonalize_ci_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L157>`_
   Replace the coefficients of the CI states by the coefficients of the
   eigenstates of the CI matrix
 
 
-`find_triples_and_quadruples <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_dress.irp.f#L206>`_
+`find_triples_and_quadruples <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_dress.irp.f#L206>`_
   Undocumented
 
 
-`h_apply_mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/H_apply.irp.f_shell_27#L416>`_
+`h_apply_mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/H_apply.irp.f_shell_27#L416>`_
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_mrcc_diexc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/H_apply.irp.f_shell_27#L1>`_
+`h_apply_mrcc_diexc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/H_apply.irp.f_shell_27#L1>`_
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_mrcc_monoexc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/H_apply.irp.f_shell_27#L268>`_
+`h_apply_mrcc_monoexc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/H_apply.irp.f_shell_27#L268>`_
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_matrix_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_utils.irp.f#L50>`_
+`h_matrix_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L50>`_
   Dressed H with Delta_ij
 
 
-`h_u_0_mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/davidson.irp.f#L360>`_
+`h_u_0_mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/davidson.irp.f#L360>`_
   Computes v_0 = H|u_0>
   .br
   n : number of determinants
@@ -122,38 +123,26 @@ Documentation
   H_jj : array of <j|H|j>
 
 
-`lambda_mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_utils.irp.f#L1>`_
+`lambda_mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L1>`_
   cm/<Psi_0|H|D_m> or perturbative 1/Delta_E(m)
 
 
-`lambda_pert <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_utils.irp.f#L2>`_
+`lambda_pert <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L2>`_
   cm/<Psi_0|H|D_m> or perturbative 1/Delta_E(m)
 
 
-`mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc.irp.f#L1>`_
+`mrcc_dress <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_dress.irp.f#L15>`_
   Undocumented
 
 
-`mrcc_dress <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_dress.irp.f#L15>`_
+`mrcc_dress_simple <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_dress.irp.f#L160>`_
   Undocumented
 
 
-`mrcc_dress_simple <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_dress.irp.f#L160>`_
+`psi_ref_lock <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_dress.irp.f#L3>`_
-  Undocumented
+  Locks on ref determinants to fill delta_ij
-
+
-
+
-`psi_cas_lock <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_dress.irp.f#L3>`_
+`run_mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_general.irp.f#L1>`_
-  Locks on CAS determinants to fill delta_ij
-
-
-`run <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc.irp.f#L10>`_
-  Undocumented
-
-
-`run_mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc.irp.f#L42>`_
-  Undocumented
-
-
-`run_mrcc_test <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc.irp.f#L29>`_
   Undocumented
 

plugins/MRCC_Utils/davidson.irp.f

@@ -35,7 +35,7 @@ subroutine davidson_diag_mrcc(dets_in,u_in,energies,dim_in,sze,N_st,Nint,iunit,i
   allocate(H_jj(sze))
   
   !$OMP PARALLEL DEFAULT(NONE)                                       &
-      !$OMP  SHARED(sze,H_jj,N_det_cas,dets_in,Nint,istate,delta_ii,idx_cas)           &
+      !$OMP  SHARED(sze,H_jj,N_det_ref,dets_in,Nint,istate,delta_ii,idx_ref)           &
       !$OMP  PRIVATE(i)
   !$OMP DO SCHEDULE(guided)
   do i=1,sze
@@ -43,8 +43,8 @@ subroutine davidson_diag_mrcc(dets_in,u_in,energies,dim_in,sze,N_st,Nint,iunit,i
   enddo
   !$OMP END DO 
   !$OMP DO SCHEDULE(guided)
-  do i=1,N_det_cas
+  do i=1,N_det_ref
-    H_jj(idx_cas(i)) +=  delta_ii(i,istate)
+    H_jj(idx_ref(i)) +=  delta_ii(i,istate)
   enddo
   !$OMP END DO 
   !$OMP END PARALLEL
@@ -384,7 +384,7 @@ subroutine H_u_0_mrcc(v_0,u_0,H_jj,n,keys_tmp,Nint,istate)
   integer, parameter             :: block_size = 157
   !$OMP PARALLEL DEFAULT(NONE)                                       &
       !$OMP PRIVATE(i,hij,j,k,idx,jj,ii,vt)                             &
-      !$OMP SHARED(n_det_cas,n_det_non_cas,idx_cas,idx_non_cas,n,H_jj,u_0,keys_tmp,Nint,v_0,istate,delta_ij)
+      !$OMP SHARED(n_det_ref,n_det_non_ref,idx_ref,idx_non_ref,n,H_jj,u_0,keys_tmp,Nint,v_0,istate,delta_ij)
   !$OMP DO SCHEDULE(static)
   do i=1,n  
     v_0(i) = H_jj(i) * u_0(i)
@@ -409,10 +409,10 @@ subroutine H_u_0_mrcc(v_0,u_0,H_jj,n,keys_tmp,Nint,istate)
   !$OMP END DO
 
   !$OMP DO SCHEDULE(guided)
-  do ii=1,n_det_cas
+  do ii=1,n_det_ref
-    i = idx_cas(ii)
+    i = idx_ref(ii)
-    do jj = 1, n_det_non_cas
+    do jj = 1, n_det_non_ref
-        j = idx_non_cas(jj)
+        j = idx_non_ref(jj)
         vt (i) = vt (i) + delta_ij(ii,jj,istate)*u_0(j)
         vt (j) = vt (j) + delta_ij(ii,jj,istate)*u_0(i)
     enddo

plugins/MRCC_Utils/mrcc_dress.irp.f

@@ -1,25 +1,25 @@
 use omp_lib
 
-BEGIN_PROVIDER [ integer(omp_lock_kind), psi_cas_lock, (psi_det_size) ]
+BEGIN_PROVIDER [ integer(omp_lock_kind), psi_ref_lock, (psi_det_size) ]
  implicit none
  BEGIN_DOC
- ! Locks on CAS determinants to fill delta_ij
+ ! Locks on ref determinants to fill delta_ij
  END_DOC
  integer :: i
  do i=1,psi_det_size
-   call omp_init_lock( psi_cas_lock(i) )
+   call omp_init_lock( psi_ref_lock(i) )
  enddo
 
 END_PROVIDER
 
-subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_cas, Ndet_non_cas,i_generator,n_selected,det_buffer,Nint,iproc)
+subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n_selected,det_buffer,Nint,iproc)
  use bitmasks
  implicit none
 
   integer, intent(in)            :: i_generator,n_selected, Nint, iproc
-  integer, intent(in) :: Ndet_cas, Ndet_non_cas
+  integer, intent(in) :: Ndet_ref, Ndet_non_ref
-  double precision, intent(inout) :: delta_ij_(Ndet_cas,Ndet_non_cas,*)
+  double precision, intent(inout) :: delta_ij_(Ndet_ref,Ndet_non_ref,*)
-  double precision, intent(inout) :: delta_ii_(Ndet_cas,*)
+  double precision, intent(inout) :: delta_ii_(Ndet_ref,*)
 
   integer(bit_kind), intent(in)  :: det_buffer(Nint,2,n_selected)
   integer                        :: i,j,k,l
@@ -43,18 +43,18 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_cas, Ndet_non_cas,i_generator,n
 
   call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq)
 
-  allocate (dIa_hla(N_states,Ndet_non_cas))
+  allocate (dIa_hla(N_states,Ndet_non_ref))
 
   ! |I>
 
   ! |alpha>
   do i_alpha=1,N_tq
-    call get_excitation_degree_vector(psi_non_cas,tq(1,1,i_alpha),degree_alpha,Nint,N_det_non_cas,idx_alpha)
+    call get_excitation_degree_vector(psi_non_ref,tq(1,1,i_alpha),degree_alpha,Nint,N_det_non_ref,idx_alpha)
 
     ! |I>
-    do i_I=1,N_det_cas
+    do i_I=1,N_det_ref
        ! Find triples and quadruple grand parents
-       call get_excitation_degree(tq(1,1,i_alpha),psi_cas(1,1,i_I),degree,Nint)
+       call get_excitation_degree(tq(1,1,i_alpha),psi_ref(1,1,i_I),degree,Nint)
        if (degree > 4) then
          cycle
        endif
@@ -65,22 +65,22 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_cas, Ndet_non_cas,i_generator,n
 
        ! <I|  <>  |alpha>
        do k_sd=1,idx_alpha(0)
-         call get_excitation_degree(psi_cas(1,1,i_I),psi_non_cas(1,1,idx_alpha(k_sd)),degree,Nint)
+         call get_excitation_degree(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),degree,Nint)
          if (degree > 2) then
            cycle
          endif
          ! <I| /k\ |alpha>
          ! <I|H|k>
-         call i_h_j(psi_cas(1,1,i_I),psi_non_cas(1,1,idx_alpha(k_sd)),Nint,hIk)
+         call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),Nint,hIk)
          do i_state=1,N_states
            dIk(i_state) = hIk * lambda_mrcc(i_state,idx_alpha(k_sd))
          enddo
          ! |l> = Exc(k -> alpha) |I>
-         call get_excitation(psi_non_cas(1,1,idx_alpha(k_sd)),tq(1,1,i_alpha),exc,degree,phase,Nint)
+         call get_excitation(psi_non_ref(1,1,idx_alpha(k_sd)),tq(1,1,i_alpha),exc,degree,phase,Nint)
          call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
          do k=1,N_int
-           tmp_det(k,1) = psi_cas(k,1,i_I)
+           tmp_det(k,1) = psi_ref(k,1,i_I)
-           tmp_det(k,2) = psi_cas(k,2,i_I)
+           tmp_det(k,2) = psi_ref(k,2,i_I)
          enddo
          ! Hole (see list_to_bitstring)
          iint = ishft(h1-1,-bit_kind_shift) + 1
@@ -108,10 +108,10 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_cas, Ndet_non_cas,i_generator,n
            dka(i_state) = 0.d0
          enddo
          do l_sd=k_sd+1,idx_alpha(0)
-           call get_excitation_degree(tmp_det,psi_non_cas(1,1,idx_alpha(l_sd)),degree,Nint)
+           call get_excitation_degree(tmp_det,psi_non_ref(1,1,idx_alpha(l_sd)),degree,Nint)
            if (degree == 0) then
-             call get_excitation(psi_cas(1,1,i_I),psi_non_cas(1,1,idx_alpha(l_sd)),exc,degree,phase2,Nint)
+             call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),exc,degree,phase2,Nint)
-             call i_h_j(psi_cas(1,1,i_I),psi_non_cas(1,1,idx_alpha(l_sd)),Nint,hIl)
+             call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hIl)
              do i_state=1,N_states
                dka(i_state) = hIl * lambda_mrcc(i_state,idx_alpha(l_sd)) * phase * phase2
              enddo
@@ -124,28 +124,28 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_cas, Ndet_non_cas,i_generator,n
        enddo
 
        do i_state=1,N_states
-         ci_inv(i_state) = 1.d0/psi_cas_coef(i_I,i_state)
+         ci_inv(i_state) = 1.d0/psi_ref_coef(i_I,i_state)
        enddo
        do l_sd=1,idx_alpha(0)
          k_sd = idx_alpha(l_sd)
-         call i_h_j(tq(1,1,i_alpha),psi_non_cas(1,1,idx_alpha(l_sd)),Nint,hla)
+         call i_h_j(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hla)
          do i_state=1,N_states
            dIa_hla(i_state,k_sd) = dIa(i_state) * hla
          enddo
        enddo
-       call omp_set_lock( psi_cas_lock(i_I) )
+       call omp_set_lock( psi_ref_lock(i_I) )
        do l_sd=1,idx_alpha(0)
          k_sd = idx_alpha(l_sd)
          do i_state=1,N_states
            delta_ij_(i_I,k_sd,i_state) += dIa_hla(i_state,k_sd)
-           if(dabs(psi_cas_coef(i_I,i_state)).ge.5.d-5)then
+           if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5)then
-            delta_ii_(i_I,i_state) -= dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_cas_coef(k_sd,i_state)
+            delta_ii_(i_I,i_state) -= dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef(k_sd,i_state)
            else
             delta_ii_(i_I,i_state)  = 0.d0
            endif
          enddo
        enddo
-       call omp_unset_lock( psi_cas_lock(i_I) )
+       call omp_unset_lock( psi_ref_lock(i_I) )
     enddo
   enddo
   deallocate (dIa_hla)
@@ -157,13 +157,13 @@ end
 
 
 
-subroutine mrcc_dress_simple(delta_ij_non_cas_,Ndet_non_cas,i_generator,n_selected,det_buffer,Nint,iproc)
+subroutine mrcc_dress_simple(delta_ij_non_ref_,Ndet_non_ref,i_generator,n_selected,det_buffer,Nint,iproc)
  use bitmasks
  implicit none
 
   integer, intent(in)            :: i_generator,n_selected, Nint, iproc
-  integer, intent(in) :: Ndet_non_cas
+  integer, intent(in) :: Ndet_non_ref
-  double precision, intent(inout) :: delta_ij_non_cas_(Ndet_non_cas,Ndet_non_cas,*)
+  double precision, intent(inout) :: delta_ij_non_ref_(Ndet_non_ref,Ndet_non_ref,*)
 
   integer(bit_kind), intent(in)  :: det_buffer(Nint,2,n_selected)
   integer                        :: i,j,k,m
@@ -184,18 +184,18 @@ subroutine mrcc_dress_simple(delta_ij_non_cas_,Ndet_non_cas,i_generator,n_select
   double precision :: f(N_states)
 
   do i=1,N_tq
-    call get_excitation_degree_vector(psi_non_cas,tq(1,1,i),degree,Nint,Ndet_non_cas,idx)
+    call get_excitation_degree_vector(psi_non_ref,tq(1,1,i),degree,Nint,Ndet_non_ref,idx)
     call i_h_j(tq(1,1,i),tq(1,1,i),Nint,haa)
     do m=1,N_states
       f(m) = 1.d0/(ci_electronic_energy(m)-haa)
     enddo
     do k=1,idx(0)
-      call i_h_j(tq(1,1,i),psi_non_cas(1,1,idx(k)),Nint,hka)
+      call i_h_j(tq(1,1,i),psi_non_ref(1,1,idx(k)),Nint,hka)
       do j=k,idx(0)
-        call i_h_j(tq(1,1,i),psi_non_cas(1,1,idx(j)),Nint,haj)
+        call i_h_j(tq(1,1,i),psi_non_ref(1,1,idx(j)),Nint,haj)
         do m=1,N_states
-          delta_ij_non_cas_(idx(k), idx(j),m) += haj*hka* f(m)
+          delta_ij_non_ref_(idx(k), idx(j),m) += haj*hka* f(m)
-          delta_ij_non_cas_(idx(j), idx(k),m) += haj*hka* f(m)
+          delta_ij_non_ref_(idx(j), idx(k),m) += haj*hka* f(m)
         enddo
       enddo 
     enddo
@@ -231,9 +231,9 @@ subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq
     endif
 
     ! Select determinants that are triple or quadruple excitations
-    ! from the CAS
+    ! from the ref
     good = .True.
-    call get_excitation_degree_vector(psi_cas,det_buffer(1,1,i),degree,Nint,N_det_cas,idx)
+    call get_excitation_degree_vector(psi_ref,det_buffer(1,1,i),degree,Nint,N_det_ref,idx)
     do k=1,idx(0)
       if (degree(k) < 3) then
         good = .False.

plugins/MRCC_Utils/mrcc_general.irp.f

@@ -0,0 +1,29 @@
+subroutine run_mrcc
+  implicit none
+  integer :: i,j
+
+  double precision :: E_new, E_old, delta_e
+  integer :: iteration
+  E_new = 0.d0
+  delta_E = 1.d0
+  iteration = 0
+  do while (delta_E > 1.d-10)
+    iteration += 1
+    print *,  '===========================' 
+    print *,  'MRCC Iteration', iteration
+    print *,  '===========================' 
+    print *,  ''
+    E_old = sum(ci_energy_dressed)
+    call write_double(6,ci_energy_dressed(1),"MRCC energy")
+    call diagonalize_ci_dressed
+    E_new = sum(ci_energy_dressed)
+    delta_E = dabs(E_new - E_old)
+    if (iteration > 20) then
+      exit
+    endif
+  enddo
+  call write_double(6,ci_energy_dressed(1),"Final MRCC energy")
+  call ezfio_set_mrcc_energy(ci_energy_dressed(1))
+  call save_wavefunction
+
+end

plugins/MRCC_Utils/mrcc_utils.irp.f

@@ -7,16 +7,16 @@
  integer :: i,k
  double precision               :: ihpsi(N_states), hii
  
- do i=1,N_det_non_cas
+ do i=1,N_det_non_ref
-   call i_h_psi(psi_non_cas(1,1,i), psi_cas, psi_cas_coef, N_int, N_det_cas,&
+   call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,&
-       size(psi_cas_coef,1), n_states, ihpsi)
+       size(psi_ref_coef,1), n_states, ihpsi)
-   call i_h_j(psi_non_cas(1,1,i),psi_non_cas(1,1,i),N_int,hii)
+   call i_h_j(psi_non_ref(1,1,i),psi_non_ref(1,1,i),N_int,hii)
    do k=1,N_states
-     lambda_pert(k,i) = 1.d0 / (psi_cas_energy_diagonalized(k)-hii)
+     lambda_pert(k,i) = 1.d0 / (psi_ref_energy_diagonalized(k)-hii)
      if (dabs(ihpsi(k)).le.1.d-3) then
        lambda_mrcc(k,i) = lambda_pert(k,i)
      else
-       lambda_mrcc(k,i) = psi_non_cas_coef(i,k)/ihpsi(k)
+       lambda_mrcc(k,i) = psi_non_ref_coef(i,k)/ihpsi(k)
      endif
    enddo
  enddo
@@ -26,17 +26,17 @@ END_PROVIDER
 
 
 
-!BEGIN_PROVIDER [ double precision, delta_ij_non_cas, (N_det_non_cas, N_det_non_cas,N_states) ]
+!BEGIN_PROVIDER [ double precision, delta_ij_non_ref, (N_det_non_ref, N_det_non_ref,N_states) ]
 !implicit none
 !BEGIN_DOC
 !! Dressing matrix in SD basis
 !END_DOC
-!delta_ij_non_cas = 0.d0
+!delta_ij_non_ref = 0.d0
-!call H_apply_mrcc_simple(delta_ij_non_cas,N_det_non_cas)
+!call H_apply_mrcc_simple(delta_ij_non_ref,N_det_non_ref)
 !END_PROVIDER
 
- BEGIN_PROVIDER [ double precision, delta_ij, (N_det_cas,N_det_non_cas,N_states) ]
+ BEGIN_PROVIDER [ double precision, delta_ij, (N_det_ref,N_det_non_ref,N_states) ]
-&BEGIN_PROVIDER [ double precision, delta_ii, (N_det_cas,N_states) ]
+&BEGIN_PROVIDER [ double precision, delta_ii, (N_det_ref,N_states) ]
  implicit none
  BEGIN_DOC
  ! Dressing matrix in N_det basis
@@ -44,7 +44,7 @@ END_PROVIDER
  integer :: i,j,m
  delta_ij = 0.d0
  delta_ii = 0.d0
- call H_apply_mrcc(delta_ij,delta_ii,N_det_cas,N_det_non_cas)
+ call H_apply_mrcc(delta_ij,delta_ii,N_det_ref,N_det_non_ref)
 END_PROVIDER
 
 BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det,N_states) ]
@@ -59,11 +59,11 @@ BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det,N_states) ]
        h_matrix_dressed(i,j,istate) = h_matrix_all_dets(i,j) 
      enddo
    enddo
-   do ii = 1, N_det_cas
+   do ii = 1, N_det_ref
-     i =idx_cas(ii)
+     i =idx_ref(ii)
      h_matrix_dressed(i,i,istate) += delta_ii(ii,istate)
-    do jj = 1, N_det_non_cas
+    do jj = 1, N_det_non_ref
-     j =idx_cas(jj)
+     j =idx_ref(jj)
      h_matrix_dressed(i,j,istate) += delta_ij(ii,jj,istate)
      h_matrix_dressed(j,i,istate) += delta_ij(ii,jj,istate)
     enddo

plugins/Psiref_CAS/.gitignore

@@ -25,5 +25,5 @@ Utils
 ezfio_interface.irp.f
 irpf90.make
 irpf90_entities
-mrcc
+mrcc_general
 tags

plugins/Psiref_CAS/NEEDED_CHILDREN_MODULES

@@ -0,0 +1 @@
+Psiref_Utils

plugins/Psiref_CAS/README.rst

@@ -0,0 +1,24 @@
+=======================
+Psiref_threshold Module
+=======================
+
+
+Reference wave function is defined as all determinants with coefficients
+such that | c_i/c_max | > threshold.
+
+Documentation
+=============
+
+.. Do not edit this section. It was auto-generated from the
+.. by the `update_README.py` script.
+
+Needed Modules
+==============
+
+.. Do not edit this section. It was auto-generated from the
+.. by the `update_README.py` script.
+
+.. image:: tree_dependency.png
+
+* `Psiref_Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Psiref_Utils>`_
+

plugins/Psiref_CAS/psi_ref.irp.f

@@ -0,0 +1,28 @@
+use bitmasks
+
+ BEGIN_PROVIDER [ integer(bit_kind), psi_ref, (N_int,2,psi_det_size) ]
+&BEGIN_PROVIDER [ double precision, psi_ref_coef,  (psi_det_size,n_states) ]
+&BEGIN_PROVIDER [ integer, idx_ref, (psi_det_size) ]
+&BEGIN_PROVIDER [ integer, N_det_ref ]
+  implicit none
+  BEGIN_DOC
+  ! CAS wave function, defined from the application of the CAS bitmask on the 
+  ! determinants. idx_cas gives the indice of the CAS determinant in psi_det.
+  END_DOC
+  integer :: i,j,k
+  N_det_ref = N_det_cas
+  do i=1,N_det_ref
+    do k=1,N_int
+      psi_ref(k,1,i) = psi_cas(k,1,i)
+      psi_ref(k,2,i) = psi_cas(k,2,i)
+    enddo
+    idx_ref(i) = idx_cas(i)
+  enddo
+  do k=1,N_states
+    do i=1,N_det_ref
+      psi_ref_coef(i,k) = psi_cas_coef(i,k)
+    enddo
+  enddo
+
+END_PROVIDER
+

plugins/Psiref_Utils/.gitignore

@@ -0,0 +1,29 @@
+# Automatically created by $QP_ROOT/scripts/module/module_handler.py 
+.ninja_deps
+.ninja_log
+AO_Basis
+Bitmask
+Determinants
+Electrons
+Ezfio_files
+Generators_full
+Hartree_Fock
+IRPF90_man
+IRPF90_temp
+Integrals_Bielec
+Integrals_Monoelec
+MOGuess
+MO_Basis
+Makefile
+Makefile.depend
+Nuclei
+Perturbation
+Properties
+Pseudo
+Selectors_full
+Utils
+ezfio_interface.irp.f
+irpf90.make
+irpf90_entities
+mrcc_general
+tags

plugins/Psiref_Utils/README.rst

@@ -0,0 +1,15 @@
+===================
+Psiref_utils Module
+===================
+
+
+Utilities related to the use of a reference wave function. This module
+needs to be loaded with any psi_ref module.
+
+
+Documentation
+=============
+
+.. Do not edit this section. It was auto-generated from the
+.. by the `update_README.py` script.
+

plugins/Psiref_Utils/psi_ref.irp.f

@@ -0,0 +1,36 @@
+use bitmasks
+
+ BEGIN_PROVIDER [ integer(bit_kind), psi_ref, (N_int,2,psi_det_size) ]
+&BEGIN_PROVIDER [ double precision, psi_ref_coef,  (psi_det_size,n_states) ]
+&BEGIN_PROVIDER [ integer, idx_ref, (psi_det_size) ]
+&BEGIN_PROVIDER [ integer, N_det_ref ]
+  implicit none
+  BEGIN_DOC
+  ! Reference wave function, defined as determinants with coefficients > 0.05
+  ! idx_ref gives the indice of the ref determinant in psi_det.
+  END_DOC
+  integer                        :: i, k, l
+  logical                        :: good
+  N_det_ref = 0
+  do i=1,N_det
+    good = .False.
+    do l = 1, N_states
+     psi_ref_coef(i,l) = 0.d0
+     good = good.or.(dabs(psi_coef(i,l)) > 0.05d0)
+    enddo
+    if (good) then
+      N_det_ref = N_det_ref+1
+      do k=1,N_int
+        psi_ref(k,1,N_det_ref) = psi_det(k,1,i)
+        psi_ref(k,2,N_det_ref) = psi_det(k,2,i)
+      enddo
+      idx_ref(N_det_ref) = i
+      do k=1,N_states
+        psi_ref_coef(N_det_ref,k) = psi_coef(i,k)
+      enddo
+    endif
+  enddo
+  call write_int(output_determinants,N_det_ref, 'Number of determinants in the reference')
+
+END_PROVIDER
+

plugins/Psiref_Utils/psi_ref_utils.irp.f

@@ -0,0 +1,123 @@
+use bitmasks
+
+
+ BEGIN_PROVIDER [ integer(bit_kind), psi_ref_sorted_bit, (N_int,2,psi_det_size) ]
+&BEGIN_PROVIDER [ double precision, psi_ref_coef_sorted_bit, (psi_det_size,N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! Reference determinants sorted to accelerate the search of a random determinant in the wave
+ ! function.
+ END_DOC
+ call sort_dets_by_det_search_key(N_det_ref, psi_ref, psi_ref_coef, &
+     psi_ref_sorted_bit, psi_ref_coef_sorted_bit)
+
+END_PROVIDER
+
+
+
+ BEGIN_PROVIDER [ integer(bit_kind), psi_non_ref,  (N_int,2,psi_det_size) ]
+&BEGIN_PROVIDER [ double precision, psi_non_ref_coef, (psi_det_size,n_states) ]
+&BEGIN_PROVIDER [ integer, idx_non_ref,  (psi_det_size) ]
+&BEGIN_PROVIDER [ integer, N_det_non_ref ]
+ implicit none
+ BEGIN_DOC
+  ! Set of determinants which are not part of the reference, defined from the application
+  ! of the reference bitmask on the determinants. 
+  ! idx_non_ref gives the indice of the determinant in psi_det.
+ END_DOC
+ integer                        :: i_non_ref,j,k
+ integer                        :: degree
+ logical                        :: in_ref
+ i_non_ref =0
+ do k=1,N_det
+   in_ref = .False.
+   do j=1,N_det_ref
+     call get_excitation_degree(psi_ref(1,1,j), psi_det(1,1,k), degree, N_int)
+     if (degree == 0) then
+       in_ref = .True.
+       exit
+     endif
+   enddo
+   if (.not.in_ref) then
+     double precision :: hij
+     i_non_ref += 1
+     do j=1,N_int
+       psi_non_ref(j,1,i_non_ref) = psi_det(j,1,k)
+       psi_non_ref(j,2,i_non_ref) = psi_det(j,2,k)
+     enddo
+     do j=1,N_states
+       psi_non_ref_coef(i_non_ref,j) = psi_coef(k,j)
+     enddo
+     idx_non_ref(i_non_ref) = k
+   endif
+ enddo
+ N_det_non_ref = i_non_ref
+END_PROVIDER
+
+ BEGIN_PROVIDER [ integer(bit_kind), psi_non_ref_sorted_bit, (N_int,2,psi_det_size) ]
+&BEGIN_PROVIDER [ double precision, psi_non_ref_coef_sorted_bit, (psi_det_size,N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! Reference determinants sorted to accelerate the search of a random determinant in the wave
+ ! function.
+ END_DOC
+ call sort_dets_by_det_search_key(N_det_ref, psi_non_ref, psi_non_ref_coef, &
+     psi_non_ref_sorted_bit, psi_non_ref_coef_sorted_bit)
+
+END_PROVIDER
+
+
+BEGIN_PROVIDER [double precision, H_matrix_ref, (N_det_ref,N_det_ref)] 
+ implicit none
+ integer :: i,j
+ double precision :: hij
+  do i = 1, N_det_ref
+   do j = 1, N_det_ref
+    call i_H_j(psi_ref(1,1,i),psi_ref(1,1,j),N_int,hij) 
+    H_matrix_ref(i,j) = hij
+   enddo
+  enddo
+END_PROVIDER
+
+ BEGIN_PROVIDER [double precision, psi_coef_ref_diagonalized, (N_det_ref,N_states)]
+&BEGIN_PROVIDER [double precision, psi_ref_energy_diagonalized, (N_states)]
+ implicit none
+ integer :: i,j
+  double precision, allocatable  :: eigenvectors(:,:), eigenvalues(:)
+  allocate (eigenvectors(size(H_matrix_ref,1),N_det_ref))
+  allocate (eigenvalues(N_det_ref))
+  call lapack_diag(eigenvalues,eigenvectors,                       &
+      H_matrix_ref,size(H_matrix_ref,1),N_det_ref)
+  do i = 1, N_states
+   psi_ref_energy_diagonalized(i) = eigenvalues(i)
+   do j = 1, N_det_ref
+    psi_coef_ref_diagonalized(j,i) = eigenvectors(j,i)
+   enddo
+  enddo
+
+
+ END_PROVIDER
+
+ BEGIN_PROVIDER [double precision, psi_ref_energy, (N_states)]
+ implicit none
+ integer :: i,j,k
+ double precision :: hij,norm,u_dot_v
+  psi_ref_energy = 0.d0
+
+
+  do k = 1, N_states
+   norm = 0.d0
+   do i = 1, N_det_ref
+    norm += psi_ref_coef(i,k) * psi_ref_coef(i,k)
+    do j = 1, N_det_ref
+      psi_ref_energy(k) += psi_ref_coef(i,k) * psi_ref_coef(j,k) * H_matrix_ref(i,j)
+    enddo
+   enddo
+   psi_ref_energy(k) = psi_ref_energy(k) /norm
+  enddo
+
+END_PROVIDER 
+
+
+
+

plugins/Psiref_threshold/.gitignore

@@ -0,0 +1,29 @@
+# Automatically created by $QP_ROOT/scripts/module/module_handler.py 
+.ninja_deps
+.ninja_log
+AO_Basis
+Bitmask
+Determinants
+Electrons
+Ezfio_files
+Generators_full
+Hartree_Fock
+IRPF90_man
+IRPF90_temp
+Integrals_Bielec
+Integrals_Monoelec
+MOGuess
+MO_Basis
+Makefile
+Makefile.depend
+Nuclei
+Perturbation
+Properties
+Pseudo
+Selectors_full
+Utils
+ezfio_interface.irp.f
+irpf90.make
+irpf90_entities
+mrcc_general
+tags

plugins/Psiref_threshold/NEEDED_CHILDREN_MODULES

@@ -0,0 +1 @@
+Psiref_Utils

plugins/Psiref_threshold/README.rst

@@ -0,0 +1,24 @@
+=======================
+Psiref_threshold Module
+=======================
+
+
+Reference wave function is defined as all determinants with coefficients
+such that | c_i/c_max | > threshold.
+
+Documentation
+=============
+
+.. Do not edit this section. It was auto-generated from the
+.. by the `update_README.py` script.
+
+Needed Modules
+==============
+
+.. Do not edit this section. It was auto-generated from the
+.. by the `update_README.py` script.
+
+.. image:: tree_dependency.png
+
+* `Psiref_Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Psiref_Utils>`_
+

plugins/Psiref_threshold/psi_ref.irp.f

@@ -0,0 +1,38 @@
+use bitmasks
+
+ BEGIN_PROVIDER [ integer(bit_kind), psi_ref, (N_int,2,psi_det_size) ]
+&BEGIN_PROVIDER [ double precision, psi_ref_coef,  (psi_det_size,n_states) ]
+&BEGIN_PROVIDER [ integer, idx_ref, (psi_det_size) ]
+&BEGIN_PROVIDER [ integer, N_det_ref ]
+  implicit none
+  BEGIN_DOC
+  ! Reference wave function, defined as determinants with coefficients > 0.05
+  ! idx_ref gives the indice of the ref determinant in psi_det.
+  END_DOC
+  integer                        :: i, k, l
+  logical                        :: good
+  double precision, parameter :: threshold=0.05d0
+  N_det_ref = 0
+  t = threshold * abs_psi_coef_max
+  do i=1,N_det
+    good = .False.
+    do l = 1, N_states
+     psi_ref_coef(i,l) = 0.d0
+     good = good.or.(dabs(psi_coef(i,l)) > t)
+    enddo
+    if (good) then
+      N_det_ref = N_det_ref+1
+      do k=1,N_int
+        psi_ref(k,1,N_det_ref) = psi_det(k,1,i)
+        psi_ref(k,2,N_det_ref) = psi_det(k,2,i)
+      enddo
+      idx_ref(N_det_ref) = i
+      do k=1,N_states
+        psi_ref_coef(N_det_ref,k) = psi_coef(i,k)
+      enddo
+    endif
+  enddo
+  call write_int(output_determinants,N_det_ref, 'Number of determinants in the reference')
+
+END_PROVIDER
+

scripts/compilation/qp_create_ninja.py

@@ -132,7 +132,7 @@ def dict_module_genelogy_path(d_module_genelogy):
 
 def get_l_module_with_ezfio_cfg():
     """
-    Return all the module who have a EZFIO.cfg
+    Return all the modules that have a EZFIO.cfg
     """
     from os import listdir
     from os.path import isfile
@@ -739,7 +739,7 @@ def create_build_ninja_module(path_module):
 
     l_string += ["rule make_all_binaries",
                  "  command = ninja -f {0}".format(ROOT_BUILD_NINJA),
-                 "  pool = console", "  description = Compile all the module",
+                 "  pool = console", "  description = Compiling all modules",
                  ""]
 
     l_string += ["rule make_clean", "  command = module_handler.py clean {0}".format(path_module.rel),
@@ -768,7 +768,7 @@ def create_build_ninja_global():
 
     l_string += ["rule make_all_binaries",
                  "  command = ninja -f {0}".format(ROOT_BUILD_NINJA),
-                 "  pool = console", "  description = Compile all the module",
+                 "  pool = console", "  description = Compiling all modules",
                  ""]
 
     l_string += ["rule make_clean",

scripts/pseudo/put_pseudo_in_ezfio.py

@@ -30,7 +30,7 @@ p = re.compile(ur'\|(\d+)><\d+\|')
 
 def get_pseudo_str(l_atom):
     """
-    Run EMSL_local for geting the str of the speudo potential
+    Run EMSL_local for getting the str of the pseudo potential
 
     str_ele :
         Element Symbol: Na

scripts/qp_set_frozen_core.py

@@ -3,7 +3,7 @@
 import os
 import sys
 
-sys.path = [ os.environ["QP_ROOT"]+"/EZFIO/Python" ] + sys.path
+sys.path = [ os.environ["QP_ROOT"]+"/install/EZFIO/Python" ] + sys.path
 from ezfio import ezfio
 
 ezfio.set_filename(sys.argv[1])

src/Determinants/README.rst

@@ -47,6 +47,14 @@ Documentation
   Needed for diag_H_mat_elem
 
 
+`abs_psi_coef_max <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L451>`_
+  Max and min values of the coefficients
+
+
+`abs_psi_coef_min <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L452>`_
+  Max and min values of the coefficients
+
+
 `ac_operator <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L1007>`_
   Needed for diag_H_mat_elem
 
@@ -120,11 +128,11 @@ Documentation
   Initial guess vectors are not necessarily orthonormal
 
 
-`connected_to_ref <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/connected_to_ref.irp.f#L155>`_
+`connected_to_ref <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/connected_to_ref.irp.f#L157>`_
   Undocumented
 
 
-`connected_to_ref_by_mono <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/connected_to_ref.irp.f#L253>`_
+`connected_to_ref_by_mono <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/connected_to_ref.irp.f#L255>`_
   Undocumented
 
 
@@ -133,7 +141,7 @@ Documentation
   After calling this subroutine, N_det, psi_det and psi_coef need to be touched
 
 
-`create_wf_of_psi_svd_matrix <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L483>`_
+`create_wf_of_psi_svd_matrix <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L416>`_
   Matrix of wf coefficients. Outer product of alpha and beta determinants
 
 
@@ -209,7 +217,7 @@ Documentation
   det_coef
 
 
-`det_connections <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L1139>`_
+`det_connections <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L1138>`_
   Build connection proxy between determinants
 
 
@@ -339,7 +347,7 @@ Documentation
   Determinants are taken from the psi_det_sorted_ab array
 
 
-`generate_all_alpha_beta_det_products <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L538>`_
+`generate_all_alpha_beta_det_products <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L471>`_
   Create a wave function from all possible alpha x beta determinants
 
 
@@ -359,15 +367,15 @@ Documentation
   Applies get_excitation_degree to an array of determinants
 
 
-`get_index_in_psi_det_alpha_unique <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L139>`_
+`get_index_in_psi_det_alpha_unique <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L141>`_
   Returns the index of the determinant in the ``psi_det_alpha_unique`` array
 
 
-`get_index_in_psi_det_beta_unique <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L221>`_
+`get_index_in_psi_det_beta_unique <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L223>`_
   Returns the index of the determinant in the ``psi_det_beta_unique`` array
 
 
-`get_index_in_psi_det_sorted_bit <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/connected_to_ref.irp.f#L48>`_
+`get_index_in_psi_det_sorted_bit <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/connected_to_ref.irp.f#L50>`_
   Returns the index of the determinant in the ``psi_det_sorted_bit`` array
 
 
@@ -383,7 +391,11 @@ Documentation
   Returns <S^2>
 
 
-`get_s2_u0 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/s2.irp.f#L82>`_
+`get_s2_u0 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/s2.irp.f#L109>`_
+  Undocumented
+
+
+`get_s2_u0_old <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/s2.irp.f#L82>`_
   Undocumented
 
 
@@ -480,7 +492,7 @@ Documentation
   .br
 
 
-`is_in_wavefunction <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/connected_to_ref.irp.f#L34>`_
+`is_in_wavefunction <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/connected_to_ref.irp.f#L36>`_
   True if the determinant ``det`` is in the wave function
 
 
@@ -500,7 +512,7 @@ Documentation
   Energy of the reference bitmask used in Slater rules
 
 
-`n_con_int <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L1131>`_
+`n_con_int <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L1130>`_
   Number of integers to represent the connections between determinants
 
 
@@ -508,11 +520,11 @@ Documentation
   Number of determinants in the wave function
 
 
-`n_det_alpha_unique <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f_template_134#L4>`_
+`n_det_alpha_unique <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f_template_136#L4>`_
   Unique alpha determinants
 
 
-`n_det_beta_unique <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f_template_134#L80>`_
+`n_det_beta_unique <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f_template_136#L80>`_
   Unique beta determinants
 
 
@@ -565,7 +577,7 @@ Documentation
   Energy of the reference bitmask used in Slater rules
 
 
-`occ_pattern_search_key <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/connected_to_ref.irp.f#L17>`_
+`occ_pattern_search_key <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/connected_to_ref.irp.f#L18>`_
   Return an integer*8 corresponding to a determinant index for searching
 
 
@@ -654,11 +666,19 @@ Documentation
   Undocumented
 
 
+`psi_coef_max <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L449>`_
+  Max and min values of the coefficients
+
+
+`psi_coef_min <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L450>`_
+  Max and min values of the coefficients
+
+
 `psi_coef_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L302>`_
   Wave function sorted by determinants contribution to the norm (state-averaged)
 
 
-`psi_coef_sorted_ab <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L450>`_
+`psi_coef_sorted_ab <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L468>`_
   Determinants on which we apply <i|H|j>.
   They are sorted by the 3 highest electrons in the alpha part,
   then by the 3 highest electrons in the beta part to accelerate
@@ -677,19 +697,19 @@ Documentation
   is empty
 
 
-`psi_det_alpha <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L25>`_
+`psi_det_alpha <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L27>`_
   List of alpha determinants of psi_det
 
 
-`psi_det_alpha_unique <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f_template_134#L3>`_
+`psi_det_alpha_unique <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f_template_136#L3>`_
   Unique alpha determinants
 
 
-`psi_det_beta <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L39>`_
+`psi_det_beta <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L41>`_
   List of beta determinants of psi_det
 
 
-`psi_det_beta_unique <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f_template_134#L79>`_
+`psi_det_beta_unique <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f_template_136#L79>`_
   Unique beta determinants
 
 
@@ -701,7 +721,7 @@ Documentation
   Wave function sorted by determinants contribution to the norm (state-averaged)
 
 
-`psi_det_sorted_ab <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L449>`_
+`psi_det_sorted_ab <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L467>`_
   Determinants on which we apply <i|H|j>.
   They are sorted by the 3 highest electrons in the alpha part,
   then by the 3 highest electrons in the beta part to accelerate
@@ -715,7 +735,7 @@ Documentation
   function.
 
 
-`psi_det_sorted_next_ab <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L451>`_
+`psi_det_sorted_next_ab <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L469>`_
   Determinants on which we apply <i|H|j>.
   They are sorted by the 3 highest electrons in the alpha part,
   then by the 3 highest electrons in the beta part to accelerate
@@ -750,31 +770,31 @@ Documentation
   psi_occ_pattern(:,2,j) = jth occ_pattern of the wave function : represent all the double occupation
 
 
-`psi_svd_alpha <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L578>`_
+`psi_svd_alpha <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L511>`_
   SVD wave function
 
 
-`psi_svd_beta <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L579>`_
+`psi_svd_beta <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L512>`_
   SVD wave function
 
 
-`psi_svd_coefs <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L580>`_
+`psi_svd_coefs <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L513>`_
   SVD wave function
 
 
-`psi_svd_matrix <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L467>`_
+`psi_svd_matrix <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L400>`_
   Matrix of wf coefficients. Outer product of alpha and beta determinants
 
 
-`psi_svd_matrix_columns <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L407>`_
+`psi_svd_matrix_columns <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L362>`_
   Matrix of wf coefficients. Outer product of alpha and beta determinants
 
 
-`psi_svd_matrix_rows <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L406>`_
+`psi_svd_matrix_rows <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L361>`_
   Matrix of wf coefficients. Outer product of alpha and beta determinants
 
 
-`psi_svd_matrix_values <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L405>`_
+`psi_svd_matrix_values <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L360>`_
   Matrix of wf coefficients. Outer product of alpha and beta determinants
 
 
@@ -782,7 +802,7 @@ Documentation
   Undocumented
 
 
-`read_dets <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L580>`_
+`read_dets <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L598>`_
   Reads the determinants from the EZFIO file
 
 
@@ -835,15 +855,15 @@ Documentation
   Save natural orbitals, obtained by diagonalization of the one-body density matrix in the MO basis
 
 
-`save_wavefunction <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L627>`_
+`save_wavefunction <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L645>`_
   Save the wave function into the EZFIO file
 
 
-`save_wavefunction_general <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L646>`_
+`save_wavefunction_general <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L664>`_
   Save the wave function into the EZFIO file
 
 
-`save_wavefunction_unsorted <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L637>`_
+`save_wavefunction_unsorted <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L655>`_
   Save the wave function into the EZFIO file
 
 
@@ -857,7 +877,7 @@ Documentation
   for a given couple of hole/particle excitations i.
 
 
-`sort_dets_by_3_highest_electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L471>`_
+`sort_dets_by_3_highest_electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L489>`_
   Determinants on which we apply <i|H|j>.
   They are sorted by the 3 highest electrons in the alpha part,
   then by the 3 highest electrons in the beta part to accelerate
@@ -890,6 +910,6 @@ Documentation
   Thresholds on selectors (fraction of the norm)
 
 
-`write_spindeterminants <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L303>`_
+`write_spindeterminants <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L305>`_
   Undocumented
 

src/Determinants/connected_to_ref.irp.f

@@ -11,6 +11,7 @@ integer*8 function det_search_key(det,Nint)
   do i=2,Nint
     det_search_key = ieor(det_search_key,iand(det(i,1),det(i,2)))
   enddo
+  det_search_key = iand(huge(det(1,1)),det_search_key)
 end
 
 
@@ -27,6 +28,7 @@ integer*8 function occ_pattern_search_key(det,Nint)
   do i=2,Nint
     occ_pattern_search_key = ieor(occ_pattern_search_key,iand(det(i,1),det(i,2)))
   enddo
+  occ_pattern_search_key = iand(huge(det(1,1)),occ_pattern_search_key)
 end
 
 

src/Determinants/determinants.irp.f

@@ -446,6 +446,24 @@ subroutine filter_3_highest_electrons( det_in, det_out, Nint )
   enddo
 end
       
+ BEGIN_PROVIDER [ double precision, psi_coef_max, (N_states) ] 
+&BEGIN_PROVIDER [ double precision, psi_coef_min, (N_states) ] 
+&BEGIN_PROVIDER [ double precision, abs_psi_coef_max, (N_states) ] 
+&BEGIN_PROVIDER [ double precision, abs_psi_coef_min, (N_states) ] 
+ implicit none
+ BEGIN_DOC
+ ! Max and min values of the coefficients
+ END_DOC
+ integer:: i
+ do i=1,N_states
+   psi_coef_min(i) = minval(psi_coef(:,i))
+   psi_coef_max(i) = maxval(psi_coef(:,i))
+   abs_psi_coef_min(i) = dabs(psi_coef_min(i))
+   abs_psi_coef_max(i) = dabs(psi_coef_max(i))
+ enddo
+
+END_PROVIDER
+
  BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_ab, (N_int,2,psi_det_size) ]
 &BEGIN_PROVIDER [ double precision, psi_coef_sorted_ab, (N_det,N_states) ]
 &BEGIN_PROVIDER [ integer, psi_det_sorted_next_ab, (2,psi_det_size) ]

src/Determinants/s2.irp.f

@@ -79,7 +79,7 @@ BEGIN_PROVIDER [ double precision, s2_values, (N_states) ]
 END_PROVIDER
 
 
-subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2)
+subroutine get_s2_u0_old(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2)
  implicit none
  use bitmasks
  integer(bit_kind), intent(in) :: psi_keys_tmp(N_int,2,nmax)
@@ -89,32 +89,60 @@ subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2)
  integer :: i,j,l
  double precision :: s2_tmp
  s2 = 0.d0
-!print*,'IN get_s2_u0'
-!print*,'n,nmax = ',n,nmax
- double precision :: accu
- accu = 0.d0
- do i = 1,n
-  accu += psi_coefs_tmp(i) * psi_coefs_tmp(i)
-! print*,'psi_coef = ',psi_coefs_tmp(i)
- enddo
-!print*,'accu = ',accu
-!print*,''
  !$OMP PARALLEL DO DEFAULT(NONE) &
  !$OMP PRIVATE(i,j,s2_tmp) SHARED(n,psi_coefs_tmp,psi_keys_tmp,N_int) REDUCTION(+:s2) SCHEDULE(dynamic) 
  do i=1,n
    do j=i+1,n
      call get_s2(psi_keys_tmp(1,1,i),psi_keys_tmp(1,1,j),s2_tmp,N_int)
-     s2 += (psi_coefs_tmp(i)+psi_coefs_tmp(i))*psi_coefs_tmp(j)*s2_tmp
+     s2 += psi_coefs_tmp(i)*psi_coefs_tmp(j)*s2_tmp
-!    s2 = s2 + 2.d0 * psi_coefs_tmp(i)*psi_coefs_tmp(j)*s2_tmp
    enddo
  enddo
  !$OMP END PARALLEL DO
+ s2 = s2+s2
  do i=1,n
    call get_s2(psi_keys_tmp(1,1,i),psi_keys_tmp(1,1,i),s2_tmp,N_int)
    s2 += psi_coefs_tmp(i)*psi_coefs_tmp(i)*s2_tmp
  enddo
  s2 +=  S_z2_Sz
-!print*,'s2 = ',s2
+end
-!print*,''
+
+subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2)
+ implicit none
+ use bitmasks
+ integer(bit_kind), intent(in) :: psi_keys_tmp(N_int,2,nmax)
+ integer, intent(in) :: n,nmax
+ double precision, intent(in) :: psi_coefs_tmp(nmax)
+ double precision, intent(out) :: s2
+ double precision :: s2_tmp
+ integer :: i,j,l,jj
+ integer, allocatable           :: idx(:)
+ s2 = 0.d0
+ !$OMP PARALLEL DEFAULT(NONE)                                        &
+     !$OMP PRIVATE(i,j,s2_tmp,idx)                                   &
+     !$OMP SHARED(n,psi_coefs_tmp,psi_keys_tmp,N_int,davidson_threshold)&
+     !$OMP REDUCTION(+:s2)
+ allocate(idx(0:n))
+ !$OMP DO SCHEDULE(dynamic)
+ do i=1,n
+   idx(0) = i
+   call filter_connected_davidson(psi_keys_tmp,psi_keys_tmp(1,1,i),N_int,i-1,idx)
+   do jj=1,idx(0)
+     j = idx(jj)
+     if ( dabs(psi_coefs_tmp(j)) + dabs(psi_coefs_tmp(i))            &
+           > davidson_threshold ) then
+       call get_s2(psi_keys_tmp(1,1,i),psi_keys_tmp(1,1,j),s2_tmp,N_int)
+       s2 = s2 + psi_coefs_tmp(i)*psi_coefs_tmp(j)*s2_tmp
+     endif
+   enddo
+ enddo
+ !$OMP END DO
+ deallocate(idx)
+ !$OMP END PARALLEL
+ s2 = s2+s2
+ do i=1,n
+   call get_s2(psi_keys_tmp(1,1,i),psi_keys_tmp(1,1,i),s2_tmp,N_int)
+   s2 = s2 + psi_coefs_tmp(i)*psi_coefs_tmp(i)*s2_tmp
+ enddo
+ s2 = s2 + S_z2_Sz
 end
 

src/Determinants/slater_rules.irp.f

@@ -1092,10 +1092,9 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
   ASSERT (Nint == N_int)
   ASSERT (n>0)
   PROVIDE ref_bitmask_energy
-  integer, parameter             :: block_size = 157
   !$OMP PARALLEL DEFAULT(NONE)                                       &
       !$OMP PRIVATE(i,hij,j,k,idx,jj,vt)                             &
-      !$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0)
+      !$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,davidson_threshold)
   !$OMP DO SCHEDULE(static)
   do i=1,n
     v_0(i) = H_jj(i) * u_0(i)
@@ -1109,7 +1108,7 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
     call filter_connected_davidson(keys_tmp,keys_tmp(1,1,i),Nint,i-1,idx)
     do jj=1,idx(0)
       j = idx(jj)
-      if ( (dabs(u_0(j)) > 1.d-7).or.((dabs(u_0(i)) > 1.d-7)) ) then
+      if ( dabs(u_0(j)) + dabs(u_0(i)) > davidson_threshold ) then
         call i_H_j(keys_tmp(1,1,j),keys_tmp(1,1,i),Nint,hij)
         vt (i) = vt (i) + hij*u_0(j)
         vt (j) = vt (j) + hij*u_0(i)

src/Determinants/spindeterminants.irp.f

@@ -14,11 +14,13 @@ integer*8 function spin_det_search_key(det,Nint)
   END_DOC
   integer, intent(in) :: Nint
   integer(bit_kind), intent(in) :: det(Nint)
+  integer(bit_kind), parameter :: unsigned_shift = not(huge(1_bit_kind)) ! 100...00
   integer :: i
   spin_det_search_key = det(1)
   do i=2,Nint
     spin_det_search_key = ieor(spin_det_search_key,det(i))
   enddo
+  spin_det_search_key = spin_det_search_key-unsigned_shift
 end
 
 
@@ -346,53 +348,6 @@ subroutine write_spindeterminants
   call ezfio_set_spindeterminants_psi_coef_matrix_rows(psi_svd_matrix_rows)
   call ezfio_set_spindeterminants_psi_coef_matrix_columns(psi_svd_matrix_columns)
   
-!  integer                        :: n_svd_coefs
-!  double precision               :: norm, f
-!  f = 1.d0/dble(N_states)
-!  norm = 1.d0
-!  do n_svd_coefs=1,N_det_alpha_unique
-!    do k=1,N_states
-!      norm -= psi_svd_coefs(n_svd_coefs,k)*psi_svd_coefs(n_svd_coefs,k)
-!    enddo
-!    if (norm < 1.d-4) then
-!      exit
-!    endif
-!  enddo
-!  n_svd_coefs -= 1
-!  call ezfio_set_spindeterminants_n_svd_coefs(n_svd_coefs)
-!  
-!  double precision, allocatable  :: dtmp(:,:,:)
-!  allocate(dtmp(N_det_alpha_unique,n_svd_coefs,N_states))
-!  do k=1,N_states
-!    do j=1,n_svd_coefs
-!      do i=1,N_det_alpha_unique
-!        dtmp(i,j,k) = psi_svd_alpha(i,j,k)
-!      enddo
-!    enddo
-!  enddo
-!  call ezfio_set_spindeterminants_psi_svd_alpha(dtmp)
-!  deallocate(dtmp)
-!
-!  allocate(dtmp(N_det_beta_unique,n_svd_coefs,N_states))
-!  do k=1,N_states
-!    do j=1,n_svd_coefs
-!      do i=1,N_det_beta_unique
-!        dtmp(i,j,k) = psi_svd_beta(i,j,k)
-!      enddo
-!    enddo
-!  enddo
-!  call ezfio_set_spindeterminants_psi_svd_beta(dtmp)
-!  deallocate(dtmp)
-!
-!  allocate(dtmp(n_svd_coefs,N_states,1))
-!  do k=1,N_states
-!    do j=1,n_svd_coefs
-!        dtmp(j,k,1) = psi_svd_coefs(j,k)
-!    enddo
-!  enddo
-!  call ezfio_set_spindeterminants_psi_svd_coefs(dtmp)
-!  deallocate(dtmp)
-
 end
 
 
@@ -419,28 +374,6 @@ BEGIN_PROVIDER  [ double precision, psi_svd_matrix_values, (N_det,N_states) ]
 
   PROVIDE psi_coef_sorted_bit
 
-!  l=0
-!  do j=1,N_det_beta_unique
-!    do k=1,N_int
-!      tmp_det(k,2) = psi_det_beta_unique(k,j)
-!    enddo
-!    do i=1,N_det_alpha_unique
-!      do k=1,N_int
-!        tmp_det(k,1) = psi_det_alpha_unique(k,i)
-!      enddo
-!      idx = get_index_in_psi_det_sorted_bit(tmp_det,N_int)
-!      if (idx > 0) then
-!        l += 1
-!        psi_svd_matrix_rows(l) = i
-!        psi_svd_matrix_columns(l) = j
-!        do k=1,N_states
-!          psi_svd_matrix_values(l,k) = psi_coef_sorted_bit(idx,k)
-!        enddo
-!      endif
-!    enddo
-!  enddo
-!  ASSERT (l == N_det)
-
   integer, allocatable :: iorder(:), to_sort(:)
   integer, external :: get_index_in_psi_det_alpha_unique
   integer, external :: get_index_in_psi_det_beta_unique

src/Integrals_Monoelec/pseudopot.f90

@@ -702,7 +702,7 @@ integer n_k(klocmax_max)
 double precision a(3),g_a,b(3),g_b,c(3),d(3)
 integer n_a(3),n_b(3),ntotA,ntotB,ntot,m
 integer i,l,k,ktot,k1,k2,k3,k1p,k2p,k3p
-double precision f,fourpi,ac,bc,freal,d2,dreal,theta_DC0,phi_DC0
+double precision f,fourpi,ac,bc,freal,d2,dreal,theta_DC0,phi_DC0,coef
 double precision,allocatable :: array_R_loc(:,:,:)
 double precision,allocatable :: array_coefs(:,:,:,:,:,:)
 double precision int_prod_bessel_loc,binom_func,accu,prod,ylm,bigI,arg
@@ -713,100 +713,119 @@ double precision int_prod_bessel_loc,binom_func,accu,prod,ylm,bigI,arg
  bc=dsqrt((b(1)-c(1))**2+(b(2)-c(2))**2+(b(3)-c(3))**2)
  arg=g_a*ac**2+g_b*bc**2
  if(arg.gt.-dlog(10.d-20))then
- Vloc=0.d0
+   Vloc=0.d0
- return
+   return
  endif
-
+ 
  ntotA=n_a(1)+n_a(2)+n_a(3)
  ntotB=n_b(1)+n_b(2)+n_b(3)
  ntot=ntotA+ntotB
-
+ 
  if(ac.eq.0.d0.and.bc.eq.0.d0)then
-  accu=0.d0
+   accu=0.d0
-
+   
-  do k=1,klocmax 
+   do k=1,klocmax
-   accu=accu+v_k(k)*crochet(n_k(k)+2+ntot,g_a+g_b+dz_k(k))
+     accu=accu+v_k(k)*crochet(n_k(k)+2+ntot,g_a+g_b+dz_k(k))
-  enddo
+   enddo
-  Vloc=accu*fourpi*bigI(0,0,0,0,n_a(1)+n_b(1),n_a(2)+n_b(2),n_a(3)+n_b(3))
+   Vloc=accu*fourpi*bigI(0,0,0,0,n_a(1)+n_b(1),n_a(2)+n_b(2),n_a(3)+n_b(3))
-  !bigI frequently is null
+   !bigI frequently is null
-  return
+   return
  endif
-
+ 
  freal=dexp(-g_a*ac**2-g_b*bc**2)
  
  d2 = 0.d0
  do i=1,3
-  d(i)=g_a*(a(i)-c(i))+g_b*(b(i)-c(i)) 
+   d(i)=g_a*(a(i)-c(i))+g_b*(b(i)-c(i))
-  d2=d2+d(i)*d(i)
+   d2=d2+d(i)*d(i)
  enddo
  d2=dsqrt(d2)
  dreal=2.d0*d2
-
+ 
-
+ 
- theta_DC0=dacos(d(3)/d2)
+ allocate (array_R_loc(-2:ntot_max+klocmax_max,klocmax_max,0:ntot_max))
- phi_DC0=datan2(d(2)/d2,d(1)/d2)
+ allocate (array_coefs(0:ntot_max,0:ntot_max,0:ntot_max,0:ntot_max,0:ntot_max,0:ntot_max))
-
+ 
- if (isnan(theta_DC0).or.isnan(phi_DC0)) then
-   print *,  'NaN in /src/Integrals_Monoelec/pseudopot.f90 at line 449.'
-   print *,  'Try to break symmetry in your molecule (1.-16 is OK).'
-   stop 1
- endif
-
-allocate (array_R_loc(-2:ntot_max+klocmax_max,klocmax_max,0:ntot_max))
-allocate (array_coefs(0:ntot_max,0:ntot_max,0:ntot_max,0:ntot_max,0:ntot_max,0:ntot_max))
-
  do ktot=-2,ntotA+ntotB+klocmax
- do l=0,ntot
+   do l=0,ntot
- do k=1,klocmax
+     do k=1,klocmax
-  array_R_loc(ktot,k,l)=freal*int_prod_bessel_loc(ktot+2,g_a+g_b+dz_k(k),l,dreal)
+       array_R_loc(ktot,k,l)=freal*int_prod_bessel_loc(ktot+2,g_a+g_b+dz_k(k),l,dreal)
- enddo
+     enddo
- enddo
- enddo
-
- do k1=0,n_a(1)
- do k2=0,n_a(2)
- do k3=0,n_a(3)
- do k1p=0,n_b(1)
- do k2p=0,n_b(2)
- do k3p=0,n_b(3)
-  array_coefs(k1,k2,k3,k1p,k2p,k3p)=binom_func(n_a(1),k1)*binom_func(n_a(2),k2)*binom_func(n_a(3),k3)  &
-  *(c(1)-a(1))**(n_a(1)-k1)*(c(2)-a(2))**(n_a(2)-k2)*(c(3)-a(3))**(n_a(3)-k3) &
-  *binom_func(n_b(1),k1p)*binom_func(n_b(2),k2p)*binom_func(n_b(3),k3p)  &
-  *(c(1)-b(1))**(n_b(1)-k1p)*(c(2)-b(2))**(n_b(2)-k2p)*(c(3)-b(3))**(n_b(3)-k3p)
- enddo
- enddo
- enddo
- enddo
- enddo
- enddo
-
- accu=0.d0
- do k=1,klocmax
-  do k1=0,n_a(1)
-  do k2=0,n_a(2)
-  do k3=0,n_a(3)
-  do k1p=0,n_b(1)
-  do k2p=0,n_b(2)
-  do k3p=0,n_b(3)
-
-  do l=0,ntot
-   do m=-l,l
-    prod=ylm(l,m,theta_DC0,phi_DC0)*array_coefs(k1,k2,k3,k1p,k2p,k3p) &
-         *bigI(l,m,0,0,k1+k1p,k2+k2p,k3+k3p)
-    ktot=k1+k2+k3+k1p+k2p+k3p+n_k(k)
-    accu=accu+prod*v_k(k)*array_R_loc(ktot,k,l)
    enddo
-  enddo
-
-  enddo
-  enddo
-  enddo
-  enddo
-  enddo
-  enddo
  enddo
+ 
+ do k1=0,n_a(1)
+   do k2=0,n_a(2)
+     do k3=0,n_a(3)
+       do k1p=0,n_b(1)
+         do k2p=0,n_b(2)
+           do k3p=0,n_b(3)
+             array_coefs(k1,k2,k3,k1p,k2p,k3p)=binom_func(n_a(1),k1)*binom_func(n_a(2),k2)*binom_func(n_a(3),k3)&
+                 *(c(1)-a(1))**(n_a(1)-k1)*(c(2)-a(2))**(n_a(2)-k2)*(c(3)-a(3))**(n_a(3)-k3)&
+                 *binom_func(n_b(1),k1p)*binom_func(n_b(2),k2p)*binom_func(n_b(3),k3p)&
+                 *(c(1)-b(1))**(n_b(1)-k1p)*(c(2)-b(2))**(n_b(2)-k2p)*(c(3)-b(3))**(n_b(3)-k3p)
+           enddo
+         enddo
+       enddo
+     enddo
+   enddo
+ enddo
+ 
+ 
+ accu=0.d0
+ if(d2 == 0.d0)then
+   l=0
+   m=0
+   coef=1.d0/dsqrt(4.d0*dacos(-1.d0))
+   do k=1,klocmax
+     do k1=0,n_a(1)
+       do k2=0,n_a(2)
+         do k3=0,n_a(3)
+           do k1p=0,n_b(1)
+             do k2p=0,n_b(2)
+               do k3p=0,n_b(3)
+                 prod=coef*array_coefs(k1,k2,k3,k1p,k2p,k3p)         &
+                     *bigI(l,m,0,0,k1+k1p,k2+k2p,k3+k3p)
+                 ktot=k1+k2+k3+k1p+k2p+k3p+n_k(k)
+                 accu=accu+prod*v_k(k)*array_R_loc(ktot,k,l)
+               enddo
+             enddo
+           enddo
+         enddo
+       enddo
+     enddo
+   enddo
+   
+ else
+   theta_DC0=dacos(d(3)/d2)
+   phi_DC0=datan2(d(2)/d2,d(1)/d2)
+   
+   do k=1,klocmax
+     do k1=0,n_a(1)
+       do k2=0,n_a(2)
+         do k3=0,n_a(3)
+           do k1p=0,n_b(1)
+             do k2p=0,n_b(2)
+               do k3p=0,n_b(3)
+                 do l=0,ntot
+                   do m=-l,l
+                     coef=ylm(l,m,theta_DC0,phi_DC0)
+                     prod=coef*array_coefs(k1,k2,k3,k1p,k2p,k3p)     &
+                         *bigI(l,m,0,0,k1+k1p,k2+k2p,k3+k3p)
+                     ktot=k1+k2+k3+k1p+k2p+k3p+n_k(k)
+                     accu=accu+prod*v_k(k)*array_R_loc(ktot,k,l)
+                   enddo
+                 enddo
+               enddo
+             enddo
+           enddo
+         enddo
+       enddo
+     enddo
+   enddo
+ endif
  Vloc=f*accu
-
+ 
  deallocate (array_R_loc)
  deallocate (array_coefs)
 end