Fix cusp

TREX-CoE · Jul 11, 2024 · 4638b40 · 4638b40
1 parent 96eea27
commit 4638b40
Show file tree

Hide file tree

Showing 2 changed files with 27 additions and 11 deletions.
diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org
@@ -5811,7 +5811,7 @@ qmckl_compute_ao_value_hpc_gaussian (const qmckl_context context,
             nucl_coord[inucl + nucl_num],
             nucl_coord[inucl + 2*nucl_num] };
 
-        /* Shift to avoid haing exact zeros at the nodes of the AOs */
+        /* Shift to avoid having exact zeros at the nodes of the AOs */
         const double shift = 1.e-20;
 
         const double x = e_coord[0] - n_coord[0];

diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org
@@ -544,23 +544,27 @@ qmckl_set_mo_basis_r_cusp(qmckl_context context,
        const double Z = qmckl_vec(ctx->nucleus.charge,inucl);
        if (Z < 0.1) continue;  // Avoid dummy atoms
        const double R = r_cusp[inucl];
+       assert (R != 0.);
+       assert (fabs(2.*R*Z-6.) > 1.e-6);
+       assert (fabs(R*R*Z-3.*R) > 1.e-6);
+       assert (fabs(R*R*R*Z-3.*R*R) > 1.e-6);
        for (int64_t i=0 ; i<ctx->mo_basis.mo_num ; ++i) {
          const double phi        = qmckl_ten3(mo_vgl_at_r_cusp_s,i,0,inucl);
          const double grad_phi   = qmckl_ten3(mo_vgl_at_r_cusp_s,i,1,inucl);
          const double lap_phi    = qmckl_ten3(mo_vgl_at_r_cusp_s,i,2,inucl);
          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.+1.e-12);
+           -(R*(2.*eta*Z-6.*grad_phi)+lap_phi*R*R+6.*phi)/(2.*R*Z-6.);
 
          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.+1.e-12);
+           (lap_phi*R*R*Z-6.*grad_phi*R*Z+6.*phi*Z+6.*eta*Z)/(2.*R*Z-6.);
 
          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+1.e-12);
+           -(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);
 
          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+1.e-12);
+           (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);
 
        }
      }
@@ -772,6 +776,7 @@ qmckl_mo_basis_select_mo (const qmckl_context context,
                           const int32_t* keep,
                           const int64_t size_max)
 {
+
   if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
     return qmckl_failwith( context,
                            QMCKL_INVALID_CONTEXT,
@@ -806,6 +811,15 @@ qmckl_mo_basis_select_mo (const qmckl_context context,
                            "Array too small: expected mo_num.");
   }
 
+/*
+  if (ctx->mo_basis.r_cusp != NULL) {
+    return qmckl_failwith( context,
+                           QMCKL_ALREADY_SET,
+                           "qmckl_get_mo_basis_select_mo",
+                           "r_cusp is already set. Please set it only after selecting MOs.");
+  }
+*/
+
   int64_t mo_num_new = 0;
   for (int64_t i=0 ; i<mo_num ; ++i) {
     if (keep[i] != 0) ++mo_num_new;
@@ -829,15 +843,17 @@ qmckl_mo_basis_select_mo (const qmckl_context context,
   ctx->mo_basis.coefficient = coefficient;
   ctx->mo_basis.mo_num = mo_num_new;
 
+  rc = qmckl_finalize_mo_basis(context);
+
   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);
+    double * tmp = (double*) calloc (ctx->nucleus.num, sizeof(double));
+    assert(tmp != NULL);
+    memcpy(tmp, ctx->mo_basis.r_cusp, ctx->nucleus.num * sizeof(double));
+    rc = qmckl_set_mo_basis_r_cusp(context, tmp, 
+        mo_num_new * 4 * ctx->nucleus.num);
+    free(tmp);
   }
 
-  rc = qmckl_finalize_mo_basis(context);
   return rc;
 }