Merge branch 'dev-lct' of https://github.com/QuantumPackage/qp2 into dev-lct

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,16 +47,47 @@ def main(arguments):
     except:
         do_pseudo = False
 
+    large = 0
+    small = 1
+
+    size = small
+    if arguments["--large"]:
+        size = large
+
     if not do_pseudo:
-        for charge in ezfio.nuclei_nucl_charge:
+
-            if charge < 5:
+        if size == large:
-                pass
+            for charge in ezfio.nuclei_nucl_charge:
-            elif charge < 13:
+		if charge <= 2:
-                n_frozen += 1
+		    pass 
-            elif charge < 31:
+		elif charge <= 10:
-                n_frozen += 5
+		    n_frozen += 1
-            else:
+		elif charge <= 18:
-                n_frozen += 9
+		    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
+                elif charge < 13:
+                    n_frozen += 1
+                elif charge < 31:
+                    n_frozen += 5
+                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))
 
 
 

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=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/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

@@ -207,61 +207,61 @@ 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"] = """
-     if (is_a_1h1p(hole).eqv..False.) cycle
+     if (.not.is_a_1h1p(hole)) cycle
     """
     self["filter_only_1h1p_double"] = """
-     if (is_a_1h1p(key).eqv..False.) cycle
+     if (.not.is_a_1h1p(key)) cycle
     """
 
 
@@ -269,22 +269,22 @@ class H_apply(object):
   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
     """
 
 
@@ -299,11 +299,11 @@ class H_apply(object):
   def filter_only_connected_to_hf(self):
     self["filter_only_connected_to_hf_single"] = """
      call connected_to_hf(hole,yes_no)
-     if (yes_no.eqv..False.) cycle
+     if (.not.yes_no) cycle
     """
     self["filter_only_connected_to_hf_double"] = """
      call connected_to_hf(key,yes_no)
-     if (yes_no.eqv..False.) cycle
+     if (.not.yes_no) cycle
     """
 
 

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
@@ -262,6 +279,100 @@ subroutine get_ao_two_e_integrals_non_zero(j,k,l,sze,out_val,out_val_index,non_z
 end
 
 
+subroutine get_ao_two_e_integrals_non_zero_jl(j,l,thresh,sze_max,sze,out_val,out_val_index,non_zero_int)
+  use map_module
+  implicit none
+  BEGIN_DOC
+  ! Gets multiple AO bi-electronic integral from the AO map .
+  ! All non-zero i are retrieved for j,k,l fixed.
+  END_DOC
+  double precision, intent(in)   :: thresh
+  integer, intent(in)            :: j,l, sze,sze_max
+  real(integral_kind), intent(out) :: out_val(sze_max)
+  integer, intent(out)           :: out_val_index(2,sze_max),non_zero_int
+
+  integer                        :: i,k
+  integer(key_kind)              :: hash
+  double precision               :: tmp
+
+  PROVIDE ao_two_e_integrals_in_map
+  non_zero_int = 0
+  if (ao_overlap_abs(j,l) < thresh) then
+    out_val = 0.d0
+    return
+  endif
+
+  non_zero_int = 0
+  do k = 1, sze
+   do i = 1, sze
+     integer, external :: ao_l4
+     double precision, external :: ao_two_e_integral
+     !DIR$ FORCEINLINE
+     if (ao_two_e_integral_schwartz(i,k)*ao_two_e_integral_schwartz(j,l) < thresh) then
+       cycle
+     endif
+     call two_e_integrals_index(i,j,k,l,hash)
+     call map_get(ao_integrals_map, hash,tmp)
+     if (dabs(tmp) < thresh ) cycle
+     non_zero_int = non_zero_int+1
+     out_val_index(1,non_zero_int) = i
+     out_val_index(2,non_zero_int) = k
+     out_val(non_zero_int) = tmp
+   enddo
+  enddo
+
+end
+
+
+subroutine get_ao_two_e_integrals_non_zero_jl_from_list(j,l,thresh,list,n_list,sze_max,out_val,out_val_index,non_zero_int)
+  use map_module
+  implicit none
+  BEGIN_DOC
+  ! Gets multiple AO two-electron integrals from the AO map .
+  ! All non-zero i are retrieved for j,k,l fixed.
+  END_DOC
+  double precision, intent(in)   :: thresh
+  integer, intent(in)            :: sze_max
+  integer, intent(in)            :: j,l, n_list,list(2,sze_max)
+  real(integral_kind), intent(out) :: out_val(sze_max)
+  integer, intent(out)           :: out_val_index(2,sze_max),non_zero_int
+
+  integer                        :: i,k
+  integer(key_kind)              :: hash
+  double precision               :: tmp
+
+  PROVIDE ao_two_e_integrals_in_map
+  non_zero_int = 0
+  if (ao_overlap_abs(j,l) < thresh) then
+    out_val = 0.d0
+    return
+  endif
+
+  non_zero_int = 0
+ integer :: kk
+  do kk = 1, n_list
+   k = list(1,kk)
+   i = list(2,kk)
+   integer, external :: ao_l4
+   double precision, external :: ao_two_e_integral
+   !DIR$ FORCEINLINE
+   if (ao_two_e_integral_schwartz(i,k)*ao_two_e_integral_schwartz(j,l) < thresh) then
+     cycle
+   endif
+   call two_e_integrals_index(i,j,k,l,hash)
+   call map_get(ao_integrals_map, hash,tmp)
+   if (dabs(tmp) < thresh ) cycle
+   non_zero_int = non_zero_int+1
+   out_val_index(1,non_zero_int) = i
+   out_val_index(2,non_zero_int) = k
+   out_val(non_zero_int) = tmp
+  enddo
+
+end
+
+
+
+
 function get_ao_map_size()
   implicit none
   integer (map_size_kind) :: get_ao_map_size

src/becke_numerical_grid/EZFIO.cfg

@@ -8,3 +8,9 @@ default: 2
 type: integer
 doc: Total number of grid points 
 interface: ezfio
+
+[thresh_grid]
+type: double precision
+doc: threshold on the weight of a given grid point
+interface: ezfio,provider,ocaml
+default: 1.e-20

src/becke_numerical_grid/grid_becke_per_atom.irp.f

@@ -0,0 +1,53 @@
+
+
+ BEGIN_PROVIDER [integer, n_pts_per_atom, (nucl_num)]
+&BEGIN_PROVIDER [integer, n_pts_max_per_atom]
+  BEGIN_DOC
+  ! Number of points which are non zero
+  END_DOC
+  integer                        :: i,j,k,l
+  n_pts_per_atom = 0
+  do j = 1, nucl_num
+    do i = 1, n_points_radial_grid -1
+      do k = 1, n_points_integration_angular
+        if(dabs(final_weight_at_r(k,i,j)) < thresh_grid)then
+          cycle
+        endif
+        n_pts_per_atom(j) += 1
+      enddo
+    enddo
+  enddo
+  n_pts_max_per_atom = maxval(n_pts_per_atom)
+END_PROVIDER
+
+ BEGIN_PROVIDER [double precision, final_grid_points_per_atom, (3,n_pts_max_per_atom,nucl_num)]
+&BEGIN_PROVIDER [double precision, final_weight_at_r_vector_per_atom, (n_pts_max_per_atom,nucl_num) ]
+&BEGIN_PROVIDER [integer, index_final_points_per_atom, (3,n_pts_max_per_atom,nucl_num) ]
+&BEGIN_PROVIDER [integer, index_final_points_per_atom_reverse, (n_points_integration_angular,n_points_radial_grid,nucl_num) ]
+ implicit none
+ integer                        :: i,j,k,l,i_count(nucl_num)
+ double precision               :: r(3)
+ i_count = 0
+  do j = 1, nucl_num
+    do i = 1, n_points_radial_grid -1
+      do k = 1, n_points_integration_angular
+        if(dabs(final_weight_at_r(k,i,j)) < thresh_grid)then
+          cycle
+        endif
+        i_count(j) += 1
+        final_grid_points_per_atom(1,i_count(j),j) = grid_points_per_atom(1,k,i,j)
+        final_grid_points_per_atom(2,i_count(j),j) = grid_points_per_atom(2,k,i,j)
+        final_grid_points_per_atom(3,i_count(j),j) = grid_points_per_atom(3,k,i,j)
+        final_weight_at_r_vector_per_atom(i_count(j),j) = final_weight_at_r(k,i,j)
+        index_final_points_per_atom(1,i_count(j),j) = k
+        index_final_points_per_atom(2,i_count(j),j) = i
+        index_final_points_per_atom(3,i_count(j),j) = j
+        index_final_points_per_atom_reverse(k,i,j) = i_count(j)
+      enddo
+    enddo
+  enddo
+
+
+
+
+END_PROVIDER 

src/becke_numerical_grid/grid_becke_vector.irp.f

@@ -8,9 +8,9 @@ BEGIN_PROVIDER [integer, n_points_final_grid]
   do j = 1, nucl_num
     do i = 1, n_points_radial_grid -1
       do k = 1, n_points_integration_angular
-!       if(dabs(final_weight_at_r(k,i,j)) < 1.d-30)then
+        if(dabs(final_weight_at_r(k,i,j)) < tresh_grid)then
-!         cycle
+          cycle
-!       endif
+        endif
         n_points_final_grid += 1
       enddo
     enddo
@@ -39,9 +39,9 @@ END_PROVIDER
   do j = 1, nucl_num
     do i = 1, n_points_radial_grid -1
       do k = 1, n_points_integration_angular
-       !if(dabs(final_weight_at_r(k,i,j)) < 1.d-30)then
+        if(dabs(final_weight_at_r(k,i,j)) < thresh_grid)then
-       !  cycle
+          cycle
-       !endif
+        endif
         i_count += 1
         final_grid_points(1,i_count) = grid_points_per_atom(1,k,i,j)
         final_grid_points(2,i_count) = grid_points_per_atom(2,k,i,j)

src/bitmask/core_inact_act_virt.irp.f

@@ -194,13 +194,13 @@ END_PROVIDER
 
 END_PROVIDER 
 
-BEGIN_PROVIDER [integer, n_inact_act ]
+BEGIN_PROVIDER [integer, n_inact_act_orb ]
  implicit none
- n_inact_act = (n_inact_orb+n_act_orb)
+ n_inact_act_orb = (n_inact_orb+n_act_orb)
 
 END_PROVIDER 
 
-BEGIN_PROVIDER [integer, list_inact_act, (n_inact_act)]
+BEGIN_PROVIDER [integer, list_inact_act, (n_inact_act_orb)]
  integer :: i,itmp
  itmp = 0
  do i = 1, n_inact_orb

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
 
 

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/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/dft_utils_in_r/ao_in_r.irp.f

@@ -121,3 +121,26 @@
  enddo
  END_PROVIDER
 
+
+ BEGIN_PROVIDER[double precision, aos_in_r_array_per_atom, (ao_num,n_pts_max_per_atom,nucl_num)]
+&BEGIN_PROVIDER[double precision, aos_in_r_array_per_atom_transp, (n_pts_max_per_atom,ao_num,nucl_num)]
+ implicit none
+ BEGIN_DOC
+ ! aos_in_r_array_per_atom(i,j,k)        = value of the ith ao on the jth grid point attached on the kth atom 
+ END_DOC
+ integer :: i,j,k
+ double precision :: aos_array(ao_num), r(3)
+ do k = 1, nucl_num
+  do i = 1, n_pts_per_atom(k)
+   r(1) = final_grid_points_per_atom(1,i,k)
+   r(2) = final_grid_points_per_atom(2,i,k)
+   r(3) = final_grid_points_per_atom(3,i,k)
+   call give_all_aos_at_r(r,aos_array)
+   do j = 1, ao_num
+    aos_in_r_array_per_atom(j,i,k) = aos_array(j)
+    aos_in_r_array_per_atom_transp(i,j,k) = aos_array(j)
+   enddo
+  enddo
+ enddo
+ END_PROVIDER
+