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

configure

@@ -213,7 +213,7 @@ EOF
               wget https://github.com/TREX-CoE/trexio/releases/download/v${VERSION}/trexio-${VERSION}.tar.gz
               tar -zxf trexio-${VERSION}.tar.gz && rm trexio-${VERSION}.tar.gz
               cd trexio-${VERSION}
-              ./configure --prefix=\${QP_ROOT} --without-hdf5
+              ./configure --prefix=\${QP_ROOT} --without-hdf5 CFLAGS='-g'
               make -j 8 && make -j 8 check && make -j 8 install
               tar -zxvf "\${QP_ROOT}"/external/qp2-dependencies/${ARCHITECTURE}/ninja.tar.gz 
               mv ninja "\${QP_ROOT}"/bin/
@@ -226,7 +226,7 @@ EOF
               wget https://github.com/TREX-CoE/trexio/releases/download/v${VERSION}/trexio-${VERSION}.tar.gz
               tar -zxf trexio-${VERSION}.tar.gz && rm trexio-${VERSION}.tar.gz
               cd trexio-${VERSION}
-              ./configure --prefix=\${QP_ROOT}
+              ./configure --prefix=\${QP_ROOT} CFLAGS="-g"
               make -j 8 && make -j 8 check && make -j 8 install
 EOF
     elif [[ ${PACKAGE} = qmckl ]] ; then
@@ -237,7 +237,7 @@ EOF
               wget https://github.com/TREX-CoE/qmckl/releases/download/v${VERSION}/qmckl-${VERSION}.tar.gz
               tar -zxf qmckl-${VERSION}.tar.gz && rm qmckl-${VERSION}.tar.gz
               cd qmckl-${VERSION}
-              ./configure --prefix=\${QP_ROOT} --enable-hpc --disable-doc
+              ./configure --prefix=\${QP_ROOT} --enable-hpc --disable-doc CFLAGS='-g'
               make && make -j 4 check && make install
 EOF
  

external/ezfio

@@ -1 +1 @@
-Subproject commit d5805497fa0ef30e70e055cde1ecec2963303e93
+Subproject commit dba01c4fe0ff7b84c5ecfb1c7c77ec68781311b3

scripts/qp_import_trexio.py

@@ -133,7 +133,7 @@ def write_ezfio(trexio_filename, filename):
     try:
         basis_type = trexio.read_basis_type(trexio_file)
 
-        if basis_type.lower() in ["gaussian", "slater"]:
+        if basis_type.lower()[0] in ["g", "s"]:
             shell_num   = trexio.read_basis_shell_num(trexio_file)
             prim_num    = trexio.read_basis_prim_num(trexio_file)
             ang_mom     = trexio.read_basis_shell_ang_mom(trexio_file)
@@ -454,19 +454,33 @@ def write_ezfio(trexio_filename, filename):
     else:
         print("None")
 
-    print("Determinant\t\t...\t", end=' ')
+    print("Determinant\t...\t", end=' ')
     alpha = [ i for i in range(num_alpha) ]
     beta  = [ i for i in range(num_beta) ]
     if trexio.has_mo_spin(trexio_file):
        spin = trexio.read_mo_spin(trexio_file)
        if max(spin) == 1:
-         beta  = [ i for i in range(mo_num) if spin[i] == 1 ]
+         alpha = [ i for i in range(len(spin)) if spin[i] == 0 ]
+         alpha = [ alpha[i] for i in range(num_alpha) ]
+         beta  = [ i for i in range(len(spin)) if spin[i] == 1 ]
          beta  = [ beta[i] for i in range(num_beta) ]
-
+         print("Warning -- UHF orbitals --", end=' ')
-    alpha = qp_bitmasks.BitMask(alpha)
+    alpha_s = ['0']*mo_num
-    beta  = qp_bitmasks.BitMask(beta )
+    beta_s  = ['0']*mo_num
-    print(alpha)
+    for i in alpha:
-    print(beta)
+      alpha_s[i] = '1'
+    for i in beta:
+      beta_s[i] = '1'
+    alpha_s = ''.join(alpha_s)[::-1]
+    beta_s = ''.join(beta_s)[::-1]
+    alpha = [ int(i,2) for i in qp_bitmasks.string_to_bitmask(alpha_s) ][::-1]
+    beta  = [ int(i,2) for i in qp_bitmasks.string_to_bitmask(beta_s ) ][::-1]
+    ezfio.set_determinants_bit_kind(8)
+    ezfio.set_determinants_n_int(1+mo_num//64)
+    ezfio.set_determinants_n_det(1)
+    ezfio.set_determinants_n_states(1)
+    ezfio.set_determinants_psi_det(alpha+beta)
+    ezfio.set_determinants_psi_coef([[1.0]])
     print("OK")
 
 

src/bi_ort_ints/three_body_ijm.irp.f

@@ -18,12 +18,13 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_direct_bi_ort, (mo_num, mo_num,
   double precision :: integral, wall1, wall0
 
   PROVIDE mo_l_coef mo_r_coef
+  provide mos_r_in_r_array_transp mos_l_in_r_array_transp
 
   three_e_3_idx_direct_bi_ort = 0.d0
   print *, ' Providing the three_e_3_idx_direct_bi_ort ...'
   call wall_time(wall0)
 
-  provide mos_r_in_r_array_transp mos_l_in_r_array_transp
+  call give_integrals_3_body_bi_ort(1, 1, 1, 1, 1, 1, integral)
 
   !$OMP PARALLEL                 &
   !$OMP DEFAULT (NONE)           &
@@ -79,6 +80,7 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_cycle_1_bi_ort, (mo_num, mo_num
 
   provide mos_r_in_r_array_transp mos_l_in_r_array_transp
 
+  call give_integrals_3_body_bi_ort(1, 1, 1, 1, 1, 1, integral)
  !$OMP PARALLEL                 &
  !$OMP DEFAULT (NONE)           &
  !$OMP PRIVATE (i,j,m,integral) & 
@@ -135,6 +137,7 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_cycle_2_bi_ort, (mo_num, mo_num
 
   provide mos_r_in_r_array_transp mos_l_in_r_array_transp
 
+  call give_integrals_3_body_bi_ort(1, 1, 1, 1, 1, 1, integral)
  !$OMP PARALLEL                 &
  !$OMP DEFAULT (NONE)           &
  !$OMP PRIVATE (i,j,m,integral) & 
@@ -191,6 +194,7 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_exch23_bi_ort, (mo_num, mo_num,
 
   provide mos_r_in_r_array_transp mos_l_in_r_array_transp
 
+  call give_integrals_3_body_bi_ort(1, 1, 1, 1, 1, 1, integral)
  !$OMP PARALLEL                 &
  !$OMP DEFAULT (NONE)           &
  !$OMP PRIVATE (i,j,m,integral) & 
@@ -247,6 +251,7 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_exch13_bi_ort, (mo_num, mo_num,
 
   provide mos_r_in_r_array_transp mos_l_in_r_array_transp
 
+  call give_integrals_3_body_bi_ort(1, 1, 1, 1, 1, 1, integral)
  !$OMP PARALLEL                 &
  !$OMP DEFAULT (NONE)           &
  !$OMP PRIVATE (i,j,m,integral) & 
@@ -303,6 +308,7 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_exch12_bi_ort, (mo_num, mo_num,
 
   provide mos_r_in_r_array_transp mos_l_in_r_array_transp
 
+  call give_integrals_3_body_bi_ort(1, 1, 1, 1, 1, 1, integral)
  !$OMP PARALLEL                 &
  !$OMP DEFAULT (NONE)           &
  !$OMP PRIVATE (i,j,m,integral) & 
@@ -349,6 +355,7 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_exch12_bi_ort_new, (mo_num, mo_
 
   provide mos_r_in_r_array_transp mos_l_in_r_array_transp
 
+  call give_integrals_3_body_bi_ort(1, 1, 1, 1, 1, 1, integral)
  !$OMP PARALLEL                 &
  !$OMP DEFAULT (NONE)           &
  !$OMP PRIVATE (i,j,m,integral) & 

src/bi_ort_ints/three_body_ints_bi_ort.irp.f

@@ -29,6 +29,9 @@ BEGIN_PROVIDER [ double precision, three_body_ints_bi_ort, (mo_num, mo_num, mo_n
 
   !provide x_W_ki_bi_ortho_erf_rk
   provide mos_r_in_r_array_transp mos_l_in_r_array_transp
+  provide int2_grad1_u12_ao_transp final_grid_points int2_grad1_u12_bimo_t
+  provide mo_l_coef mo_r_coef mos_l_in_r_array_transp mos_r_in_r_array_transp n_points_final_grid
+
 
  !$OMP PARALLEL                       &
  !$OMP DEFAULT (NONE)                 &

src/non_h_ints_mu/jast_deriv.irp.f

@@ -121,76 +121,93 @@
     double precision :: f
     f = 1.d0 / dble(elec_num - 1)
 
+    integer*8 :: n_points, k
+    n_points = n_points_extra_final_grid * n_points_final_grid
+
     double precision, allocatable :: rij(:,:,:)
-    allocate( rij(3, 2, n_points_extra_final_grid) )
+    allocate( rij(3, 2, n_points) )
 
     use qmckl
     integer(qmckl_exit_code) :: rc
 
-    integer*8 :: npoints
-    npoints = n_points_extra_final_grid
-
     double precision, allocatable :: gl(:,:,:)
-    allocate( gl(2,4,n_points_extra_final_grid) )
+    allocate( gl(2,4,n_points) )
 
+    k=0
     do ipoint = 1, n_points_final_grid  ! r1
-
       do jpoint = 1, n_points_extra_final_grid ! r2
-        rij(1:3, 1, jpoint) = final_grid_points      (1:3, ipoint)
+        k=k+1
-        rij(1:3, 2, jpoint) = final_grid_points_extra(1:3, jpoint)
+        rij(1:3, 1, k) = final_grid_points      (1:3, ipoint)
+        rij(1:3, 2, k) = final_grid_points_extra(1:3, jpoint)
       enddo
+    enddo
 
 
-      rc = qmckl_set_electron_coord(qmckl_ctx_jastrow, 'N', npoints, rij, npoints*6_8)
+    rc = qmckl_set_electron_coord(qmckl_ctx_jastrow, 'N', n_points, rij, n_points*6_8)
-      if (rc /= QMCKL_SUCCESS) then
+    if (rc /= QMCKL_SUCCESS) then
-        print *, irp_here, 'qmckl error in set_electron_coord'
+      print *, irp_here, 'qmckl error in set_electron_coord'
-        stop -1
+      rc = qmckl_check(qmckl_ctx_jastrow, rc)
-      endif
+      stop -1
+    endif
 
 
-      ! ---
+    ! ---
-      ! e-e term
+    ! e-e term
 
-      rc = qmckl_get_jastrow_champ_factor_ee_gl(qmckl_ctx_jastrow, gl, 8_8*npoints)
+    rc = qmckl_get_jastrow_champ_factor_ee_gl(qmckl_ctx_jastrow, gl, 8_8*n_points)
-      if (rc /= QMCKL_SUCCESS) then
+    if (rc /= QMCKL_SUCCESS) then
-        print *, irp_here, 'qmckl error in fact_ee_gl'
+      print *, irp_here, ' qmckl error in fact_ee_gl'
-        stop -1
+      rc = qmckl_check(qmckl_ctx_jastrow, rc)
-      endif
+      stop -1
+    endif
 
+    k=0
+    do ipoint = 1, n_points_final_grid  ! r1
       do jpoint = 1, n_points_extra_final_grid ! r2
-        grad1_u12_num(jpoint,ipoint,1) = gl(1,1,jpoint)
+        k=k+1
-        grad1_u12_num(jpoint,ipoint,2) = gl(1,2,jpoint)
+        grad1_u12_num(jpoint,ipoint,1) = gl(1,1,k)
-        grad1_u12_num(jpoint,ipoint,3) = gl(1,3,jpoint)
+        grad1_u12_num(jpoint,ipoint,2) = gl(1,2,k)
+        grad1_u12_num(jpoint,ipoint,3) = gl(1,3,k)
       enddo
+    enddo
 
     ! ---
     ! e-e-n term
 
-!    rc = qmckl_get_jastrow_champ_factor_een_gl(qmckl_ctx_jastrow, gl, 8_8*npoints)
+!    rc = qmckl_get_jastrow_champ_factor_een_gl(qmckl_ctx_jastrow, gl, 8_8*n_points)
 !    if (rc /= QMCKL_SUCCESS) then
 !      print *, irp_here, 'qmckl error in fact_een_gl'
+!        rc = qmckl_check(qmckl_ctx_jastrow, rc)
 !      stop -1
 !    endif
 !
+!    k=0
+!    do ipoint = 1, n_points_final_grid  ! r1
 !    do jpoint = 1, n_points_extra_final_grid ! r2
-!        grad1_u12_num(jpoint,ipoint,1) = grad1_u12_num(jpoint,ipoint,1) + gl(1,1,jpoint)
+!        k=k+1
-!        grad1_u12_num(jpoint,ipoint,2) = grad1_u12_num(jpoint,ipoint,2) + gl(1,2,jpoint)
+!        grad1_u12_num(jpoint,ipoint,1) = grad1_u12_num(jpoint,ipoint,1) + gl(1,1,k)
-!        grad1_u12_num(jpoint,ipoint,3) = grad1_u12_num(jpoint,ipoint,3) + gl(1,3,jpoint)
+!        grad1_u12_num(jpoint,ipoint,2) = grad1_u12_num(jpoint,ipoint,2) + gl(1,2,k)
+!        grad1_u12_num(jpoint,ipoint,3) = grad1_u12_num(jpoint,ipoint,3) + gl(1,3,k)
+!    enddo
 !    enddo
 
       ! ---
       ! e-n term
 
-      rc = qmckl_get_jastrow_champ_factor_en_gl(qmckl_ctx_jastrow, gl, 8_8*npoints)
+    rc = qmckl_get_jastrow_champ_factor_en_gl(qmckl_ctx_jastrow, gl, 8_8*n_points)
-      if (rc /= QMCKL_SUCCESS) then
+    if (rc /= QMCKL_SUCCESS) then
-        print *, irp_here, 'qmckl error in fact_en_gl'
+      print *, irp_here, 'qmckl error in fact_en_gl'
-        stop -1
+      rc = qmckl_check(qmckl_ctx_jastrow, rc)
-      endif
+      stop -1
+    endif
 
+    k=0
+    do ipoint = 1, n_points_final_grid  ! r1
       do jpoint = 1, n_points_extra_final_grid ! r2
-        grad1_u12_num(jpoint,ipoint,1) = grad1_u12_num(jpoint,ipoint,1) + f * gl(1,1,jpoint)
+        k = k+1
-        grad1_u12_num(jpoint,ipoint,2) = grad1_u12_num(jpoint,ipoint,2) + f * gl(1,2,jpoint)
+        grad1_u12_num(jpoint,ipoint,1) = grad1_u12_num(jpoint,ipoint,1) + f * gl(1,1,k)
-        grad1_u12_num(jpoint,ipoint,3) = grad1_u12_num(jpoint,ipoint,3) + f * gl(1,3,jpoint)
+        grad1_u12_num(jpoint,ipoint,2) = grad1_u12_num(jpoint,ipoint,2) + f * gl(1,2,k)
+        grad1_u12_num(jpoint,ipoint,3) = grad1_u12_num(jpoint,ipoint,3) + f * gl(1,3,k)
       enddo
 
       do jpoint = 1, n_points_extra_final_grid ! r2
@@ -199,7 +216,6 @@
         dz = grad1_u12_num(jpoint,ipoint,3)
         grad1_u12_squared_num(jpoint,ipoint) = dx*dx + dy*dy + dz*dz
       enddo
-
     enddo
 
     deallocate(gl, rij)

src/non_h_ints_mu/qmckl.irp.f

@@ -59,6 +59,10 @@ BEGIN_PROVIDER [ integer*8, qmckl_ctx_jastrow ]
   if (rc /= QMCKL_SUCCESS) stop -1
 
 
+  rc =  qmckl_set_jastrow_champ_cord_num(qmckl_ctx_jastrow, jast_qmckl_cord_num*1_8)
+  rc = qmckl_check(qmckl_ctx_jastrow, rc)
+  if (rc /= QMCKL_SUCCESS) stop -1
+
   if (jast_qmckl_cord_num > 0) then
     rc =  qmckl_set_jastrow_champ_c_vector(qmckl_ctx_jastrow, jast_qmckl_c_vector, 1_8*jast_qmckl_c_vector_size)
     rc = qmckl_check(qmckl_ctx_jastrow, rc)

src/tc_bi_ortho/tc_bi_ortho.irp.f

@@ -13,6 +13,17 @@ program tc_bi_ortho
   my_n_pt_a_grid = tc_grid1_a
   touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
 
+  if(j1b_type .ge. 100) then
+    my_extra_grid_becke  = .True.
+    PROVIDE tc_grid2_a tc_grid2_r
+    my_n_pt_r_extra_grid = tc_grid2_r
+    my_n_pt_a_extra_grid = tc_grid2_a
+    touch my_extra_grid_becke my_n_pt_r_extra_grid my_n_pt_a_extra_grid
+
+    call write_int(6, my_n_pt_r_extra_grid, 'radial  internal grid over')
+    call write_int(6, my_n_pt_a_extra_grid, 'angular internal grid over')
+  endif
+
   read_wf = .True.
   touch read_wf
 

src/tc_scf/tc_petermann_factor.irp.f

@@ -30,9 +30,22 @@ subroutine main()
 
   allocate(Sl(mo_num,mo_num), Sr(mo_num,mo_num), Pf(mo_num,mo_num))
 
-  call dgemm( "T", "N", mo_num, mo_num, ao_num, 1.d0                       &
+
-            , mo_l_coef, size(mo_l_coef, 1), mo_l_coef, size(mo_l_coef, 1) &
+  call LTxSxR(ao_num, mo_num, mo_l_coef, ao_overlap, mo_r_coef, Sl)
-            , 0.d0, Sl, size(Sl, 1) )
+  !call dgemm( "T", "N", mo_num, mo_num, ao_num, 1.d0                       &
+  !          , mo_l_coef, size(mo_l_coef, 1), mo_l_coef, size(mo_l_coef, 1) &
+  !          , 0.d0, Sl, size(Sl, 1) )
+
+  print *, ''
+  print *, ' left-right orthog matrix:'
+  do i = 1, mo_num
+    write(*,'(100(F8.4,X))') Sl(:,i)
+  enddo
+
+  call LTxSxR(ao_num, mo_num, mo_l_coef, ao_overlap, mo_l_coef, Sl)
+  !call dgemm( "T", "N", mo_num, mo_num, ao_num, 1.d0                       &
+  !          , mo_l_coef, size(mo_l_coef, 1), mo_l_coef, size(mo_l_coef, 1) &
+  !          , 0.d0, Sl, size(Sl, 1) )
 
   print *, ''
   print *, ' left-orthog matrix:'
@@ -40,9 +53,10 @@ subroutine main()
     write(*,'(100(F8.4,X))') Sl(:,i)
   enddo
 
-  call dgemm( "T", "N", mo_num, mo_num, ao_num, 1.d0                       &
+  call LTxSxR(ao_num, mo_num, mo_r_coef, ao_overlap, mo_r_coef, Sr)
-            , mo_r_coef, size(mo_r_coef, 1), mo_r_coef, size(mo_r_coef, 1) &
+!  call dgemm( "T", "N", mo_num, mo_num, ao_num, 1.d0                       &
-            , 0.d0, Sr, size(Sr, 1) )
+!            , mo_r_coef, size(mo_r_coef, 1), mo_r_coef, size(mo_r_coef, 1) &
+!            , 0.d0, Sr, size(Sr, 1) )
 
   print *, ''
   print *, ' right-orthog matrix:'