Skip to content

Commit

Permalink
Fixed bug in cusp
Browse files Browse the repository at this point in the history
  • Loading branch information
@scemama
scemama committed Jul 11, 2024
1 parent 9553059 commit 96eea27
Show file tree
Hide file tree
Showing 2 changed files with 65 additions and 59 deletions.
4 changes: 2 additions & 2 deletions org/qmckl_blas.org
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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");
}

Expand Down
120 changes: 63 additions & 57 deletions org/qmckl_mo.org
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);

}
}
Expand Down Expand Up @@ -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;
}
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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

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
! 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

end do
end if
end do ! inucl
end do
info = QMCKL_SUCCESS

Expand Down Expand Up @@ -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];
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) )));
Expand Down

0 comments on commit 96eea27

Please sign in to comment.