TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-12-22 12:23:56 +01:00
Code Issues Projects Releases Wiki Activity

Fix cusp

Browse Source
This commit is contained in:
Anthony Scemama 2024-07-11 13:53:10 +02:00
parent 96eea27713
commit 4638b40724
2 changed files with 27 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
org/qmckl_ao.org
View File

@ -5811,7 +5811,7 @@ qmckl_compute_ao_value_hpc_gaussian (const qmckl_context context,
nucl_coord[inucl + nucl_num], nucl_coord[inucl + nucl_num],
nucl_coord[inucl + 2*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 shift = 1.e-20;
const double x = e_coord[0] - n_coord[0]; const double x = e_coord[0] - n_coord[0];

36
org/qmckl_mo.org
View File

@ -544,6 +544,10 @@ qmckl_set_mo_basis_r_cusp(qmckl_context context,
const double Z = qmckl_vec(ctx->nucleus.charge,inucl); const double Z = qmckl_vec(ctx->nucleus.charge,inucl);
if (Z < 0.1) continue; // Avoid dummy atoms if (Z < 0.1) continue; // Avoid dummy atoms
const double R = r_cusp[inucl]; 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) { 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 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 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); const double eta = qmckl_mat(mo_value_at_nucl_no_s,i,inucl);
qmckl_ten3(ctx->mo_basis.cusp_param,i,0,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) = 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) = 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) = 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 int32_t* keep,
const int64_t size_max) const int64_t size_max)
{ {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( context, return qmckl_failwith( context,
QMCKL_INVALID_CONTEXT, QMCKL_INVALID_CONTEXT,
@ -806,6 +811,15 @@ qmckl_mo_basis_select_mo (const qmckl_context context,
"Array too small: expected mo_num."); "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; int64_t mo_num_new = 0;
for (int64_t i=0 ; i<mo_num ; ++i) { for (int64_t i=0 ; i<mo_num ; ++i) {
if (keep[i] != 0) ++mo_num_new; 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.coefficient = coefficient;
ctx->mo_basis.mo_num = mo_num_new; ctx->mo_basis.mo_num = mo_num_new;
rc = qmckl_finalize_mo_basis(context);
if (ctx->mo_basis.r_cusp != NULL) { if (ctx->mo_basis.r_cusp != NULL) {
double * r_cusp_old = calloc(ctx->nucleus.num, sizeof(double)); double * tmp = (double*) calloc (ctx->nucleus.num, sizeof(double));
assert (r_cusp_old != NULL); assert(tmp != NULL);
memcpy(r_cusp_old, ctx->mo_basis.r_cusp, ctx->nucleus.num*sizeof(double)); memcpy(tmp, ctx->mo_basis.r_cusp, ctx->nucleus.num * sizeof(double));
qmckl_set_mo_basis_r_cusp(context, r_cusp_old, ctx->nucleus.num); rc = qmckl_set_mo_basis_r_cusp(context, tmp,
free(r_cusp_old); mo_num_new * 4 * ctx->nucleus.num);
free(tmp);
} }
rc = qmckl_finalize_mo_basis(context);
return rc; return rc;
} }