1
0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2025-04-30 04:15:00 +02:00
This commit is contained in:
Emiel Slootman 2025-02-07 09:30:09 +01:00
parent 45efa5cb5e
commit 88e4a6b48c
4 changed files with 99 additions and 85 deletions

View File

@ -319,10 +319,13 @@ typedef struct qmckl_ao_basis_struct {
uint64_t primitive_vgl_date;
double * restrict shell_vgl;
uint64_t shell_vgl_date;
uint64_t shell_vgl_maxsize;
double * restrict ao_vgl;
uint64_t ao_vgl_date;
uint64_t ao_vgl_maxsize;
double * restrict ao_value;
uint64_t ao_value_date;
uint64_t ao_value_maxsize;
double * restrict shell_hessian;
uint64_t shell_hessian_date;
@ -4051,17 +4054,11 @@ qmckl_exit_code qmckl_provide_ao_basis_shell_vgl(qmckl_context context)
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->ao_basis.shell_num * 5 * ctx->point.num * sizeof(double);
if (ctx->ao_basis.shell_vgl != NULL) {
qmckl_memory_info_struct mem_info_test = qmckl_memory_info_struct_zero;
rc = qmckl_get_malloc_info(context, ctx->ao_basis.shell_vgl, &mem_info_test);
/* if rc != QMCKL_SUCCESS, we are maybe in an _inplace function because the
memory was not allocated with qmckl_malloc */
if ((rc == QMCKL_SUCCESS) && (mem_info_test.size != mem_info.size)) {
rc = qmckl_free(context, ctx->ao_basis.shell_vgl);
assert (rc == QMCKL_SUCCESS);
ctx->ao_basis.shell_vgl = NULL;
if (mem_info.size > ctx->ao_basis.shell_vgl_maxsize) {
if (ctx->ao_basis.shell_vgl != NULL) {
rc = qmckl_free(context, ctx->ao_basis.shell_vgl);
assert (rc == QMCKL_SUCCESS);
ctx->ao_basis.shell_vgl = NULL;
}
}
@ -4069,7 +4066,7 @@ qmckl_exit_code qmckl_provide_ao_basis_shell_vgl(qmckl_context context)
if (ctx->ao_basis.shell_vgl == NULL) {
double* shell_vgl = (double*) qmckl_malloc(context, mem_info);
ctx->ao_basis.shell_vgl_maxsize = mem_info.size;
if (shell_vgl == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
@ -6742,24 +6739,19 @@ qmckl_exit_code qmckl_provide_ao_basis_ao_value(qmckl_context context)
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->ao_basis.ao_num * ctx->point.num * sizeof(double);
if (ctx->ao_basis.ao_value != NULL) {
qmckl_memory_info_struct mem_info_test = qmckl_memory_info_struct_zero;
rc = qmckl_get_malloc_info(context, ctx->ao_basis.ao_value, &mem_info_test);
/* if rc != QMCKL_SUCCESS, we are maybe in an _inplace function because the
memory was not allocated with qmckl_malloc */
if ((rc == QMCKL_SUCCESS) && (mem_info_test.size != mem_info.size)) {
rc = qmckl_free(context, ctx->ao_basis.ao_value);
assert (rc == QMCKL_SUCCESS);
ctx->ao_basis.ao_value = NULL;
}
if (mem_info.size > ctx->ao_basis.ao_value_maxsize) {
if (ctx->ao_basis.ao_value != NULL) {
rc = qmckl_free(context, ctx->ao_basis.ao_value);
assert (rc == QMCKL_SUCCESS);
ctx->ao_basis.ao_value = NULL;
}
}
/* Allocate array */
if (ctx->ao_basis.ao_value == NULL) {
double* ao_value = (double*) qmckl_malloc(context, mem_info);
ctx->ao_basis.ao_value_maxsize = mem_info.size;
if (ao_value == NULL) {
return qmckl_failwith( context,
@ -7650,25 +7642,20 @@ qmckl_exit_code qmckl_provide_ao_basis_ao_vgl(qmckl_context context)
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->ao_basis.ao_num * 5 * ctx->point.num * sizeof(double);
if (ctx->ao_basis.ao_vgl != NULL) {
qmckl_memory_info_struct mem_info_test = qmckl_memory_info_struct_zero;
rc = qmckl_get_malloc_info(context, ctx->ao_basis.ao_vgl, &mem_info_test);
/* if rc != QMCKL_SUCCESS, we are maybe in an _inplace function because the
memory was not allocated with qmckl_malloc */
if ((rc == QMCKL_SUCCESS) && (mem_info_test.size != mem_info.size)) {
rc = qmckl_free(context, ctx->ao_basis.ao_vgl);
assert (rc == QMCKL_SUCCESS);
ctx->ao_basis.ao_vgl = NULL;
}
if (mem_info.size > ctx->ao_basis.ao_vgl_maxsize) {
if (ctx->ao_basis.ao_vgl != NULL) {
rc = qmckl_free(context, ctx->ao_basis.ao_vgl);
assert (rc == QMCKL_SUCCESS);
ctx->ao_basis.ao_vgl = NULL;
}
}
/* Allocate array */
if (ctx->ao_basis.ao_vgl == NULL) {
double* ao_vgl = (double*) qmckl_malloc(context, mem_info);
ctx->ao_basis.ao_vgl_maxsize = mem_info.size;
if (ao_vgl == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,

View File

@ -100,8 +100,10 @@ typedef struct qmckl_forces_struct{
uint64_t forces_jastrow_een_l_date;
double * restrict forces_ao_value;
uint64_t forces_ao_value_date;
uint64_t forces_ao_value_maxsize;
double * restrict forces_mo_value;
uint64_t forces_mo_value_date;
uint64_t forces_mo_value_maxsize;
double * restrict forces_mo_g;
uint64_t forces_mo_g_date;
double * restrict forces_mo_l;
@ -4362,27 +4364,27 @@ integer(qmckl_exit_code) function qmckl_compute_forces_jastrow_delta_p_doc( &
do nw=1, walk_num
do m=0, cord_num-1
do m=0, cord_num-1
do j = 1, elec_num
delta_e(j) = een_rescaled_single_e(j,m,nw) - een_rescaled_e(j,num,m,nw)
end do
do j = 1, elec_num
delta_e(j) = een_rescaled_single_e(j,m,nw) - een_rescaled_e(j,num,m,nw)
end do
do l=0, cord_num
do k = 1, 3
do a = 1, nucl_num
delta_n_gl = een_rescaled_single_n_gl(k,a,l,nw) - een_rescaled_n_gl(num, k, a, l, nw)
do j = 1, elec_num
forces_delta_p(j,k,a,l,m,nw) = delta_e(j) * (-een_rescaled_single_n_gl(k,a,l,nw)) - &
een_rescaled_e(j,num,m,nw) * delta_n_gl
end do
do j = 1, elec_num
forces_delta_p(num,k,a,l,m,nw) = forces_delta_p(num,k,a,l,m,nw) - delta_e(j) * een_rescaled_n_gl(j,k,a,l,nw)
end do
end do
end do
do l=0, cord_num
do a = 1, nucl_num
do k = 1, 3
delta_n_gl = een_rescaled_single_n_gl(k,a,l,nw) - een_rescaled_n_gl(num, k, a, l, nw)
do j = 1, elec_num
forces_delta_p(j,k,a,l,m,nw) = delta_e(j) * (-een_rescaled_single_n_gl(k,a,l,nw)) - &
een_rescaled_e(j,num,m,nw) * delta_n_gl
end do
do j = 1, elec_num
forces_delta_p(num,k,a,l,m,nw) = forces_delta_p(num,k,a,l,m,nw) - delta_e(j) * een_rescaled_n_gl(j,k,a,l,nw)
end do
end do
end do
end do
end do
end do
end do
end function qmckl_compute_forces_jastrow_delta_p_doc
@ -5308,24 +5310,23 @@ qmckl_exit_code qmckl_provide_forces_ao_value(qmckl_context context)
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->ao_basis.ao_num * 3 * ctx->nucleus.num * ctx->point.num * sizeof(double);
if (ctx->forces.forces_ao_value != NULL) {
qmckl_memory_info_struct mem_info_test = qmckl_memory_info_struct_zero;
rc = qmckl_get_malloc_info(context, ctx->forces.forces_ao_value, &mem_info_test);
/* if rc != QMCKL_SUCCESS, we are maybe in an _inplace function because the
memory was not allocated with qmckl_malloc */
if ((rc == QMCKL_SUCCESS) && (mem_info_test.size != mem_info.size)) {
rc = qmckl_free(context, ctx->forces.forces_ao_value);
assert (rc == QMCKL_SUCCESS);
ctx->forces.forces_ao_value = NULL;
if (mem_info.size > ctx->forces.forces_ao_value_maxsize) {
if (ctx->forces.forces_ao_value != NULL) {
rc = qmckl_free(context, ctx->forces.forces_ao_value);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context, rc,
"qmckl_provide_forces_ao_value",
"Unable to free ctx->forces.forces_ao_value");
}
ctx->forces.forces_ao_value = NULL;
}
}
/* Allocate array */
if (ctx->forces.forces_ao_value == NULL) {
double* forces_ao_value = (double*) qmckl_malloc(context, mem_info);
ctx->forces.forces_ao_value_maxsize = mem_info.size;
if (forces_ao_value == NULL) {
return qmckl_failwith( context,
@ -5670,14 +5671,8 @@ qmckl_exit_code qmckl_provide_forces_mo_value(qmckl_context context)
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = 3 * ctx->mo_basis.mo_num * ctx->point.num * ctx->nucleus.num * sizeof(double);
if (ctx->forces.forces_mo_value != NULL) {
qmckl_memory_info_struct mem_info_test = qmckl_memory_info_struct_zero;
rc = qmckl_get_malloc_info(context, ctx->forces.forces_mo_value, &mem_info_test);
/* if rc != QMCKL_SUCCESS, we are maybe in an _inplace function because the
memory was not allocated with qmckl_malloc */
if ((rc == QMCKL_SUCCESS) && (mem_info_test.size != mem_info.size)) {
if (mem_info.size > ctx->forces.forces_mo_value_maxsize) {
if (ctx->forces.forces_mo_value != NULL) {
rc = qmckl_free(context, ctx->forces.forces_mo_value);
assert (rc == QMCKL_SUCCESS);
ctx->forces.forces_mo_value = NULL;

View File

View File

@ -717,6 +717,7 @@ qmckl_exit_code qmckl_compute_single_ee_distance (
#+begin_src c :tangle (eval c_test)
printf("Single e-e distance\n");
/* Check if Jastrow is properly initialized */
assert(qmckl_jastrow_champ_provided(context));
@ -764,6 +765,8 @@ for (int elec = 0; elec < elec_num; elec++){
}
}
}
printf("OK\n");
#+end_src
** Electron-nucleus distances
@ -979,6 +982,8 @@ qmckl_exit_code qmckl_compute_single_en_distance (
#+begin_src c :tangle (eval c_test)
printf("Single e-n distance\n");
/* Check if Jastrow is properly initialized */
assert(qmckl_jastrow_champ_provided(context));
@ -1020,6 +1025,8 @@ for (int elec = 0; elec < elec_num; elec++){
}
}
}
printf("OK\n");
#+end_src
* Electron-electron-nucleus Jastrow
@ -1377,6 +1384,8 @@ qmckl_compute_een_rescaled_single_e (const qmckl_context context,
#+begin_src c :tangle (eval c_test)
printf("Single een rescaled e-e distance\n");
/* Check if Jastrow is properly initialized */
assert(qmckl_jastrow_champ_provided(context));
@ -1424,6 +1433,8 @@ for (int elec = 0; elec < elec_num; elec++){
}
}
printf("OK\n");
#+end_src
@ -1666,6 +1677,7 @@ end function qmckl_compute_een_rescaled_single_n
#+begin_src c :tangle (eval c_test)
printf("Single een rescaled e-n distance\n");
/* Check if Jastrow is properly initialized */
assert(qmckl_jastrow_champ_provided(context));
@ -1711,6 +1723,8 @@ for (int elec = 0; elec < elec_num; elec++){
}
}
}
printf("OK\n");
#+end_src
** $\delta P$ matrix
@ -1982,6 +1996,7 @@ qmckl_compute_jastrow_champ_delta_p (const qmckl_context context,
#+begin_src c :tangle (eval c_test)
printf("Delta p\n");
/* Check if Jastrow is properly initialized */
assert(qmckl_jastrow_champ_provided(context));
@ -2032,6 +2047,8 @@ for (int elec = 0; elec < elec_num; elec++){
}
}
}
printf("OK\n");
#+end_src
** Electron-electron-nucleus Jastrow value
@ -2354,7 +2371,8 @@ qmckl_compute_jastrow_champ_factor_single_een (const qmckl_context context,
*** Test :noexport:
#+begin_src c :tangle (eval c_test)
printf("Delta een\n");
/* Check if Jastrow is properly initialized */
assert(qmckl_jastrow_champ_provided(context));
@ -2402,6 +2420,8 @@ for (int elec = 0; elec < elec_num; elec++){
}
printf("OK\n");
#+end_src
** Electron-nucleus rescaled distance derivative
@ -2920,9 +2940,6 @@ integer function qmckl_compute_een_rescaled_single_e_gl_doc( &
elec_dist_gl(4, i) = 2.0d0 * rij_inv
end do
! print *, single_ee_distance(2, 1), elec_dist_gl(1,2)
! stop
elec_dist_gl(:, num) = 0.0d0
do l = 1, cord_num
@ -3668,7 +3685,7 @@ integer(qmckl_exit_code) function qmckl_compute_jastrow_champ_factor_single_een_
double precision :: een_rescaled_delta_n_gl(4, nucl_num, 0:cord_num, walk_num)
double precision :: een_rescaled_delta_n(nucl_num, 0:cord_num, walk_num)
double precision :: dpg1_m, dpg1_ml, dp_m, dp_ml, een_r_m, een_r_ml, een_r_gl_m, een_r_gl_ml
num = num_in + 1
info = QMCKL_SUCCESS
@ -3699,13 +3716,28 @@ integer(qmckl_exit_code) function qmckl_compute_jastrow_champ_factor_single_een_
do a = 1, nucl_num
cn = c_vector_full(a, n)
if(cn == 0.d0) cycle
!do i = 1, elec_num
! delta_een_gl(i,kk,nw) = delta_een_gl(i,kk,nw) + ( &
! delta_p_gl(i,a,kk,m ,k,nw) * een_rescaled_n(i,a,m+l,nw) + &
! delta_p_gl(i,a,kk,m+l,k,nw) * een_rescaled_n(i,a,m ,nw) + &
! delta_p(i,a,m ,k,nw) * een_rescaled_n_gl(i,kk,a,m+l,nw) + &
! delta_p(i,a,m+l,k,nw) * een_rescaled_n_gl(i,kk,a,m ,nw) ) * cn
!end do
do i = 1, elec_num
delta_een_gl(i,kk,nw) = delta_een_gl(i,kk,nw) + ( &
delta_p_gl(i,a,kk,m ,k,nw) * een_rescaled_n(i,a,m+l,nw) + &
delta_p_gl(i,a,kk,m+l,k,nw) * een_rescaled_n(i,a,m ,nw) + &
delta_p(i,a,m ,k,nw) * een_rescaled_n_gl(i,kk,a,m+l,nw) + &
delta_p(i,a,m+l,k,nw) * een_rescaled_n_gl(i,kk,a,m ,nw) ) * cn
end do
! Cache repeated accesses
dpg1_m = delta_p_gl(i,a,kk,m ,k,nw)
dpg1_ml = delta_p_gl(i,a,kk,m+l,k,nw)
dp_m = delta_p(i,a,m ,k,nw)
dp_ml = delta_p(i,a,m+l,k,nw)
een_r_m = een_rescaled_n(i,a,m ,nw)
een_r_ml = een_rescaled_n(i,a,m+l,nw)
een_r_gl_m = een_rescaled_n_gl(i,kk,a,m ,nw)
een_r_gl_ml = een_rescaled_n_gl(i,kk,a,m+l,nw)
delta_een_gl(i,kk,nw) = delta_een_gl(i,kk,nw) + cn * &
(dpg1_m * een_r_ml + dpg1_ml * een_r_m + dp_m * een_r_gl_ml + dp_ml * een_r_gl_m)
end do
delta_een_gl(num,kk,nw) = delta_een_gl(num,kk,nw) + ( &
(dtmp_c(num,kk,a,m ,k,nw) + delta_p_gl(num,a,kk,m ,k,nw)) * een_rescaled_delta_n(a,m+l,nw) + &