From 4638b4072428c949b8787723b9a8f1b23f352be3 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 11 Jul 2024 13:53:10 +0200 Subject: [PATCH] Fix cusp --- org/qmckl_ao.org | 2 +- org/qmckl_mo.org | 36 ++++++++++++++++++++++++++---------- 2 files changed, 27 insertions(+), 11 deletions(-) diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org index 29322ae1..284daa1d 100644 --- 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 index 3339c929..ed456ad3 100644 --- a/org/qmckl_mo.org +++ b/org/qmckl_mo.org @@ -544,6 +544,10 @@ 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 ; imo_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); @@ -551,16 +555,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.+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 ; imo_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; }