From 96eea277137c4f4b36c97e877790533d0b2b4f5e Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 11 Jul 2024 10:54:21 +0200 Subject: [PATCH] Fixed bug in cusp --- org/qmckl_blas.org | 4 +- org/qmckl_mo.org | 118 ++++++++++++++++++++++++--------------------- 2 files changed, 64 insertions(+), 58 deletions(-) diff --git a/org/qmckl_blas.org b/org/qmckl_blas.org index c9a4545..596f941 100644 --- a/org/qmckl_blas.org +++ b/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"); } diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org index 33f5432..3339c92 100644 --- a/org/qmckl_mo.org +++ b/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,51 +2671,55 @@ 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 - inucl = ao_nucl(k)+1 - if ( (en_distance(inucl,j) < r_cusp(inucl)) .and. (ao_ang_mom(k) == 0) ) cycle - c1 = ao_vgl(k,1,j) - c2 = ao_vgl(k,2,j) - c3 = ao_vgl(k,3,j) - c4 = ao_vgl(k,4,j) - c5 = ao_vgl(k,5,j) - do i=1,mo_num - mo_vgl(i,1,j) = mo_vgl(i,1,j) + coefficient_t(i,k) * c1 - mo_vgl(i,2,j) = mo_vgl(i,2,j) + coefficient_t(i,k) * c2 - mo_vgl(i,3,j) = mo_vgl(i,3,j) + coefficient_t(i,k) * c3 - 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 + if (ao_vgl(k,1,j) /= 0.d0) then + inucl = ao_nucl(k)+1 + 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) + c4 = ao_vgl(k,4,j) + c5 = ao_vgl(k,5,j) + do i=1,mo_num + mo_vgl(i,1,j) = mo_vgl(i,1,j) + coefficient_t(i,k) * c1 + mo_vgl(i,2,j) = mo_vgl(i,2,j) + coefficient_t(i,k) * c2 + mo_vgl(i,3,j) = mo_vgl(i,3,j) + coefficient_t(i,k) * c3 + 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 + ! Cusp adjustment + do inucl=1,nucl_num + r = en_distance(inucl,j) + 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,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) + + mo_vgl(i,2,j) = mo_vgl(i,2,j) + r_vec(1) * c1 + mo_vgl(i,3,j) = mo_vgl(i,3,j) + r_vec(2) * c1 + mo_vgl(i,4,j) = mo_vgl(i,4,j) + r_vec(3) * c1 + + mo_vgl(i,5,j) = mo_vgl(i,5,j) + & + 2.d0*cusp_param(i,2,inucl)*r_inv + & + 6.d0*cusp_param(i,3,inucl) + & + 12.d0*cusp_param(i,4,inucl)*r - 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,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) - - mo_vgl(i,2,j) = mo_vgl(i,2,j) + r_vec(1) * c1 - mo_vgl(i,3,j) = mo_vgl(i,3,j) + r_vec(2) * c1 - mo_vgl(i,4,j) = mo_vgl(i,4,j) + r_vec(3) * c1 - - mo_vgl(i,5,j) = mo_vgl(i,5,j) + & - 2.d0*cusp_param(i,2,inucl)*r_inv + & - 6.d0*cusp_param(i,3,inucl) + & - 12.d0*cusp_param(i,4,inucl)*r - - enddo - enddo ! inucl + end do + 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