Fixed bug in cusp

org/qmckl_blas.org

@@ -719,14 +719,14 @@ qmckl_vector_of_double(const qmckl_context context,
   if (vector.size == 0) {
     return qmckl_failwith( context,
                            QMCKL_INVALID_ARG_4,
-                           "qmckl_double_of_vector",
+                           "qmckl_vector_of_double",
                            "Vector not allocated");
   }
 
   if (vector.size != size_max) {
       return qmckl_failwith( context,
                              QMCKL_INVALID_ARG_4,
-                             "qmckl_double_of_vector",
+                             "qmckl_vector_of_double",
                              "Wrong vector size");
   }
 

org/qmckl_mo.org

@@ -551,16 +551,16 @@ qmckl_set_mo_basis_r_cusp(qmckl_context context,
          const double eta        = qmckl_mat(mo_value_at_nucl_no_s,i,inucl);
 
          qmckl_ten3(ctx->mo_basis.cusp_param,i,0,inucl) =
-           -(R*(2.*eta*Z-6.*grad_phi)+lap_phi*R*R+6.*phi)/(2.*R*Z-6.);
+           -(R*(2.*eta*Z-6.*grad_phi)+lap_phi*R*R+6.*phi)/(2.*R*Z-6.+1.e-12);
 
          qmckl_ten3(ctx->mo_basis.cusp_param,i,1,inucl) =
-           (lap_phi*R*R*Z-6.*grad_phi*R*Z+6.*phi*Z+6.*eta*Z)/(2.*R*Z-6.);
+           (lap_phi*R*R*Z-6.*grad_phi*R*Z+6.*phi*Z+6.*eta*Z)/(2.*R*Z-6.+1.e-12);
 
          qmckl_ten3(ctx->mo_basis.cusp_param,i,2,inucl) =
-           -(R*(-5.*grad_phi*Z-1.5*lap_phi)+lap_phi*R*R*Z+3.*phi*Z+3.*eta*Z+6.*grad_phi)/(R*R*Z-3.*R);
+           -(R*(-5.*grad_phi*Z-1.5*lap_phi)+lap_phi*R*R*Z+3.*phi*Z+3.*eta*Z+6.*grad_phi)/(R*R*Z-3.*R+1.e-12);
 
          qmckl_ten3(ctx->mo_basis.cusp_param,i,3,inucl) =
-           (R*(-2.*grad_phi*Z-lap_phi)+0.5*lap_phi*R*R*Z+phi*Z+eta*Z+3.*grad_phi)/(R*R*R*Z-3.*R*R);
+           (R*(-2.*grad_phi*Z-lap_phi)+0.5*lap_phi*R*R*Z+phi*Z+eta*Z+3.*grad_phi)/(R*R*R*Z-3.*R*R+1.e-12);
 
        }
      }
@@ -829,6 +829,14 @@ qmckl_mo_basis_select_mo (const qmckl_context context,
   ctx->mo_basis.coefficient = coefficient;
   ctx->mo_basis.mo_num = mo_num_new;
 
+  if (ctx->mo_basis.r_cusp != NULL) {
+    double * r_cusp_old = calloc(ctx->nucleus.num, sizeof(double));
+    assert (r_cusp_old != NULL);
+    memcpy(r_cusp_old, ctx->mo_basis.r_cusp, ctx->nucleus.num*sizeof(double));
+    qmckl_set_mo_basis_r_cusp(context, r_cusp_old, ctx->nucleus.num);
+    free(r_cusp_old);
+  }
+
   rc = qmckl_finalize_mo_basis(context);
   return rc;
 }
@@ -1825,12 +1833,6 @@ integer function qmckl_compute_mo_basis_mo_vgl_doc_f(context, &
      end do
   end do
 
-! print *, 'coucou'
-! info = qmckl_dgemm(context,'N', 'N', mo_num, point_num*5_8, ao_num, 1.d0, &
-!      coefficient_t, int(size(coefficient_t,1),8),      &
-!      ao_vgl, int(size(ao_vgl,1),8), 0.d0,                  &
-!      mo_vgl, int(size(mo_vgl,1),8))
-
 end function qmckl_compute_mo_basis_mo_vgl_doc_f
     #+end_src
 
@@ -2669,9 +2671,9 @@ integer function qmckl_compute_mo_basis_mo_vgl_cusp_doc_f(context, &
      ! Initial contribution of the MO
      mo_vgl(:,:,j) = 0.d0
      do k=1,ao_num
-        if (ao_vgl(k,1,j) == 0.d0) cycle
+        if (ao_vgl(k,1,j) /= 0.d0) then
           inucl = ao_nucl(k)+1
-        if ( (en_distance(inucl,j) < r_cusp(inucl)) .and. (ao_ang_mom(k) == 0) ) cycle
+          if ( (en_distance(inucl,j) > r_cusp(inucl)) .or. (ao_ang_mom(k) > 0) ) then
              c1 = ao_vgl(k,1,j)
              c2 = ao_vgl(k,2,j)
              c3 = ao_vgl(k,3,j)
@@ -2684,21 +2686,24 @@ integer function qmckl_compute_mo_basis_mo_vgl_cusp_doc_f(context, &
                 mo_vgl(i,4,j) = mo_vgl(i,4,j) + coefficient_t(i,k) * c4
                 mo_vgl(i,5,j) = mo_vgl(i,5,j) + coefficient_t(i,k) * c5
              end do
+          end if
+       end if
      end do
 
     ! Cusp adjustment
     do inucl=1,nucl_num
        r = en_distance(inucl,j)
-        if (r > r_cusp(inucl)) cycle
+       if (r < r_cusp(inucl)) then
           
           r_vec(1:3) = point_coord(j,1:3) - nucl_coord(inucl,1:3)
           r_inv = 1.d0/r
           
-
           do i=1,mo_num
              mo_vgl(i,1,j) = mo_vgl(i,1,j) +  &
-                cusp_param(i,1,inucl) + r*(cusp_param(i,2,inucl) + r*(  &
+                  cusp_param(i,1,inucl) + r*( &
-                cusp_param(i,3,inucl) + r* cusp_param(i,4,inucl)     ))
+                  cusp_param(i,2,inucl) + r*( &
+                  cusp_param(i,3,inucl) + r*( &
+                  cusp_param(i,4,inucl)    )))
              
              c1 = r_inv * cusp_param(i,2,inucl) + 2.d0*cusp_param(i,3,inucl) +  &
                   r * 3.d0 * cusp_param(i,4,inucl)
@@ -2712,8 +2717,9 @@ integer function qmckl_compute_mo_basis_mo_vgl_cusp_doc_f(context, &
                   6.d0*cusp_param(i,3,inucl) +       &
                   12.d0*cusp_param(i,4,inucl)*r
 
-        enddo
+          end do
-     enddo ! inucl
+       end if
+    end do ! inucl
   end do
   info = QMCKL_SUCCESS
 
@@ -3029,17 +3035,17 @@ qmckl_compute_mo_basis_mo_vgl_cusp_hpc (const qmckl_context context,
 
     // TODO
     for (int64_t inucl=0 ; inucl<nucl_num ; ++inucl) {
-      if (ria[inucl] < r_cusp[inucl]) {
       const double r = ria[inucl];
+      if (r < r_cusp[inucl]) {
         const double r_vec[3] = {
           qmckl_mat(point_coord,ipoint,0) - qmckl_mat(nucl_coord,inucl,0),
           qmckl_mat(point_coord,ipoint,1) - qmckl_mat(nucl_coord,inucl,1),
           qmckl_mat(point_coord,ipoint,2) - qmckl_mat(nucl_coord,inucl,2) };
         const double r_inv = 1./r;
 
-IVDEP
         for (int64_t i=0 ; i<mo_num ; ++i) {
-          vgl1[i] = vgl1[i] + qmckl_ten3(cusp_param,i,0,inucl) + r*(
+          vgl1[i] = vgl1[i] +
+                     qmckl_ten3(cusp_param,i,0,inucl) + r*(
                      qmckl_ten3(cusp_param,i,1,inucl) + r*(
                      qmckl_ten3(cusp_param,i,2,inucl) + r*(
                      qmckl_ten3(cusp_param,i,3,inucl) )));