mirror of
https://github.com/TREX-CoE/qmckl.git
synced 2024-07-18 08:53:47 +02:00
Finished een_rescale_n. #22
This commit is contained in:
vijay gopal chilkuri
parent
a7c1fe526e
commit
b6bb9be359
@ -2993,10 +2993,12 @@ integer function qmckl_compute_een_rescaled_e_f(context, walk_num, elec_num, cor
|
||||
return
|
||||
endif
|
||||
|
||||
! Prepare table of exponentiated distances raised to appropriate power
|
||||
een_rescaled_e = 0.0d0
|
||||
do nw = 1, walk_num
|
||||
een_rescaled_e_ij = 0.0d0
|
||||
een_rescaled_e_ij(:, 1) = 1.0d0
|
||||
|
||||
! Prepare table of exponentiated distances raised to appropriate power
|
||||
do nw = 1, walk_num
|
||||
k = 0
|
||||
do j = 1, elec_num
|
||||
do i = 1, j - 1
|
||||
@ -3007,13 +3009,12 @@ integer function qmckl_compute_een_rescaled_e_f(context, walk_num, elec_num, cor
|
||||
|
||||
do l = 2, cord_num
|
||||
do k = 1, elec_num * (elec_num - 1)/2
|
||||
een_rescaled_e_ij(k, l + 1) = een_rescaled_e_ij(k, l + 1 - 1) * een_rescaled_e_ij(k, 1)
|
||||
een_rescaled_e_ij(k, l + 1) = een_rescaled_e_ij(k, l + 1 - 1) * een_rescaled_e_ij(k, 2)
|
||||
end do
|
||||
end do
|
||||
|
||||
! prepare the actual een table
|
||||
een_rescaled_e = 0.0d0
|
||||
een_rescaled_e(0, :, :, :) = 1.0d0
|
||||
een_rescaled_e(0, :, :, nw) = 1.0d0
|
||||
do l = 1, cord_num
|
||||
k = 0
|
||||
do j = 1, elec_num
|
||||
@ -3090,7 +3091,7 @@ een_rescaled_e_ij[:,0] = 1.0
|
||||
|
||||
k = 0
|
||||
for j in range(elec_num):
|
||||
for i in range(j - 1):
|
||||
for i in range(j):
|
||||
een_rescaled_e_ij[k, 1] = np.exp(-kappa * elec_dist[i, j])
|
||||
k = k + 1
|
||||
|
||||
@ -3104,7 +3105,7 @@ een_rescaled_e[:,:,0] = 1.0
|
||||
for l in range(1,cord_num+1):
|
||||
k = 0
|
||||
for j in range(elec_num):
|
||||
for i in range(j - 1):
|
||||
for i in range(j):
|
||||
x = een_rescaled_e_ij[k, l]
|
||||
een_rescaled_e[i, j, l] = x
|
||||
een_rescaled_e[j, i, l] = x
|
||||
@ -3134,18 +3135,314 @@ double een_rescaled_e[walk_num][elec_num][elec_num][(cord_num + 1)];
|
||||
rc = qmckl_get_jastrow_een_rescaled_e(context, &(een_rescaled_e[0][0][0][0]));
|
||||
|
||||
// value of (0,2,1)
|
||||
//printf("%10.15f = \n", een_rescaled_e[0][0][2][1]);
|
||||
//printf("%10.15f = \n", een_rescaled_e[0][0][3][1]);
|
||||
//printf("%10.15f = \n", een_rescaled_e[0][0][4][1]);
|
||||
//printf("%10.15f = \n", een_rescaled_e[0][1][3][2]);
|
||||
//printf("%10.15f = \n", een_rescaled_e[0][1][4][2]);
|
||||
//printf("%10.15f = \n", een_rescaled_e[0][1][5][2]);
|
||||
//assert(fabs(een_rescaled_e[0][0][2][1]-) < 1.e-12);
|
||||
//assert(fabs(een_rescaled_e[0][0][3][1]-) < 1.e-12);
|
||||
//assert(fabs(een_rescaled_e[0][0][4][1]-) < 1.e-12);
|
||||
//assert(fabs(een_rescaled_e[0][1][3][2]-) < 1.e-12);
|
||||
//assert(fabs(een_rescaled_e[0][1][4][2]-) < 1.e-12);
|
||||
//assert(fabs(een_rescaled_e[0][1][5][2]-) < 1.e-12);
|
||||
assert(fabs(een_rescaled_e[0][0][2][1]-0.08084493981483197) < 1.e-12);
|
||||
assert(fabs(een_rescaled_e[0][0][3][1]-0.1066745707571846) < 1.e-12);
|
||||
assert(fabs(een_rescaled_e[0][0][4][1]-0.01754273169464735) < 1.e-12);
|
||||
assert(fabs(een_rescaled_e[0][1][3][2]-0.02214680362033448) < 1.e-12);
|
||||
assert(fabs(een_rescaled_e[0][1][4][2]-0.0005700154999202759) < 1.e-12);
|
||||
assert(fabs(een_rescaled_e[0][1][5][2]-0.3424402276009091) < 1.e-12);
|
||||
|
||||
#+end_src
|
||||
|
||||
** Electron-nucleus rescaled distances for each order
|
||||
|
||||
~een_rescaled_n~ stores the table of the rescaled distances between
|
||||
electrons and nucleii raised to the power \(p\) defined by ~cord_num~:
|
||||
|
||||
\[
|
||||
C_{ia,p} = \left( 1 - \exp{-\kappa C_{ia}} \right)^p
|
||||
\]
|
||||
|
||||
where \(C_{ia}\) is the matrix of electron-nucleus distances.
|
||||
|
||||
*** Get
|
||||
|
||||
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
|
||||
qmckl_exit_code qmckl_get_jastrow_een_rescaled_n(qmckl_context context, double* const distance_rescaled);
|
||||
#+end_src
|
||||
|
||||
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
||||
qmckl_exit_code qmckl_get_jastrow_een_rescaled_n(qmckl_context context, double* const distance_rescaled)
|
||||
{
|
||||
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
||||
return QMCKL_NULL_CONTEXT;
|
||||
}
|
||||
|
||||
qmckl_exit_code rc;
|
||||
|
||||
rc = qmckl_provide_een_rescaled_n(context);
|
||||
if (rc != QMCKL_SUCCESS) return rc;
|
||||
|
||||
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
||||
assert (ctx != NULL);
|
||||
|
||||
size_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1);
|
||||
memcpy(distance_rescaled, ctx->jastrow.een_rescaled_n, sze * sizeof(double));
|
||||
|
||||
return QMCKL_SUCCESS;
|
||||
}
|
||||
#+end_src
|
||||
|
||||
*** Provide :noexport:
|
||||
|
||||
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
|
||||
qmckl_exit_code qmckl_provide_een_rescaled_n(qmckl_context context);
|
||||
#+end_src
|
||||
|
||||
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
||||
qmckl_exit_code qmckl_provide_een_rescaled_n(qmckl_context context)
|
||||
{
|
||||
|
||||
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
||||
return QMCKL_NULL_CONTEXT;
|
||||
}
|
||||
|
||||
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
||||
assert (ctx != NULL);
|
||||
|
||||
/* Check if ee distance is provided */
|
||||
qmckl_exit_code rc = qmckl_provide_en_distance(context);
|
||||
if(rc != QMCKL_SUCCESS) return rc;
|
||||
|
||||
/* Compute if necessary */
|
||||
if (ctx->date > ctx->jastrow.een_rescaled_n_date) {
|
||||
|
||||
/* Allocate array */
|
||||
if (ctx->jastrow.een_rescaled_n == NULL) {
|
||||
|
||||
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
||||
mem_info.size = ctx->electron.num * ctx->nucleus.num *
|
||||
ctx->electron.walk_num * (ctx->jastrow.cord_num + 1) * sizeof(double);
|
||||
double* een_rescaled_n = (double*) qmckl_malloc(context, mem_info);
|
||||
|
||||
if (een_rescaled_n == NULL) {
|
||||
return qmckl_failwith( context,
|
||||
QMCKL_ALLOCATION_FAILED,
|
||||
"qmckl_een_rescaled_n",
|
||||
NULL);
|
||||
}
|
||||
ctx->jastrow.een_rescaled_n = een_rescaled_n;
|
||||
}
|
||||
|
||||
qmckl_exit_code rc =
|
||||
qmckl_compute_een_rescaled_n(context,
|
||||
ctx->electron.walk_num,
|
||||
ctx->electron.num,
|
||||
ctx->nucleus.num,
|
||||
ctx->jastrow.cord_num,
|
||||
ctx->electron.rescale_factor_kappa_en,
|
||||
ctx->electron.en_distance,
|
||||
ctx->jastrow.een_rescaled_n);
|
||||
if (rc != QMCKL_SUCCESS) {
|
||||
return rc;
|
||||
}
|
||||
|
||||
ctx->jastrow.een_rescaled_n_date = ctx->date;
|
||||
}
|
||||
|
||||
return QMCKL_SUCCESS;
|
||||
}
|
||||
#+end_src
|
||||
|
||||
*** Compute
|
||||
:PROPERTIES:
|
||||
:Name: qmckl_compute_een_rescaled_n
|
||||
:CRetType: qmckl_exit_code
|
||||
:FRetType: qmckl_exit_code
|
||||
:END:
|
||||
|
||||
#+NAME: qmckl_factor_een_rescaled_n_args
|
||||
| qmckl_context | context | in | Global state |
|
||||
| int64_t | walk_num | in | Number of walkers |
|
||||
| int64_t | elec_num | in | Number of electrons |
|
||||
| int64_t | nucl_num | in | Number of atoms |
|
||||
| int64_t | cord_num | in | Order of polynomials |
|
||||
| double | rescale_factor_kappa_en | in | Factor to rescale ee distances |
|
||||
| double | en_distance[walk_num][elec_num][nucl_num] | in | Electron-nucleus distances |
|
||||
| double | een_rescaled_n[walk_num][elec_num][nucl_num][0:cord_num] | out | Electron-nucleus rescaled distances |
|
||||
|
||||
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
||||
integer function qmckl_compute_een_rescaled_n_f(context, walk_num, elec_num, nucl_num, cord_num, rescale_factor_kappa_en, &
|
||||
en_distance, een_rescaled_n) &
|
||||
result(info)
|
||||
use qmckl
|
||||
implicit none
|
||||
integer(qmckl_context), intent(in) :: context
|
||||
integer*8 , intent(in) :: walk_num
|
||||
integer*8 , intent(in) :: elec_num
|
||||
integer*8 , intent(in) :: nucl_num
|
||||
integer*8 , intent(in) :: cord_num
|
||||
double precision , intent(in) :: rescale_factor_kappa_en
|
||||
double precision , intent(in) :: en_distance(elec_num,nucl_num,walk_num)
|
||||
double precision , intent(out) :: een_rescaled_n(0:cord_num,nucl_num,elec_num,walk_num)
|
||||
double precision :: x
|
||||
integer*8 :: i, a, k, l, nw
|
||||
|
||||
info = QMCKL_SUCCESS
|
||||
|
||||
if (context == QMCKL_NULL_CONTEXT) then
|
||||
info = QMCKL_INVALID_CONTEXT
|
||||
return
|
||||
endif
|
||||
|
||||
if (walk_num <= 0) then
|
||||
info = QMCKL_INVALID_ARG_2
|
||||
return
|
||||
endif
|
||||
|
||||
if (elec_num <= 0) then
|
||||
info = QMCKL_INVALID_ARG_3
|
||||
return
|
||||
endif
|
||||
|
||||
if (nucl_num <= 0) then
|
||||
info = QMCKL_INVALID_ARG_4
|
||||
return
|
||||
endif
|
||||
|
||||
if (cord_num <= 0) then
|
||||
info = QMCKL_INVALID_ARG_5
|
||||
return
|
||||
endif
|
||||
|
||||
! Prepare table of exponentiated distances raised to appropriate power
|
||||
een_rescaled_n = 0.0d0
|
||||
do nw = 1, walk_num
|
||||
|
||||
! prepare the actual een table
|
||||
een_rescaled_n(0, :, :, nw) = 1.0d0
|
||||
|
||||
do a = 1, nucl_num
|
||||
do i = 1, elec_num
|
||||
een_rescaled_n(1, a, i, nw) = dexp(-rescale_factor_kappa_en * en_distance(i, a, nw))
|
||||
end do
|
||||
end do
|
||||
|
||||
do l = 2, cord_num
|
||||
do a = 1, nucl_num
|
||||
do i = 1, elec_num
|
||||
een_rescaled_n(l, a, i, nw) = een_rescaled_n(l - 1, a, i, nw) * een_rescaled_n(1, a, i, nw)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
end function qmckl_compute_een_rescaled_n_f
|
||||
#+end_src
|
||||
|
||||
#+CALL: generate_c_header(table=qmckl_factor_een_rescaled_n_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
||||
|
||||
#+RESULTS:
|
||||
#+begin_src c :tangle (eval h_func) :comments org
|
||||
qmckl_exit_code qmckl_compute_een_rescaled_n (
|
||||
const qmckl_context context,
|
||||
const int64_t walk_num,
|
||||
const int64_t elec_num,
|
||||
const int64_t nucl_num,
|
||||
const int64_t cord_num,
|
||||
const double rescale_factor_kappa_en,
|
||||
const double* en_distance,
|
||||
double* const een_rescaled_n );
|
||||
#+end_src
|
||||
|
||||
#+CALL: generate_c_interface(table=qmckl_factor_een_rescaled_n_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
|
||||
|
||||
#+RESULTS:
|
||||
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
||||
integer(c_int32_t) function qmckl_compute_een_rescaled_n &
|
||||
(context, &
|
||||
walk_num, &
|
||||
elec_num, &
|
||||
nucl_num, &
|
||||
cord_num, &
|
||||
rescale_factor_kappa_en, &
|
||||
en_distance, &
|
||||
een_rescaled_n) &
|
||||
bind(C) result(info)
|
||||
|
||||
use, intrinsic :: iso_c_binding
|
||||
implicit none
|
||||
|
||||
integer (c_int64_t) , intent(in) , value :: context
|
||||
integer (c_int64_t) , intent(in) , value :: walk_num
|
||||
integer (c_int64_t) , intent(in) , value :: elec_num
|
||||
integer (c_int64_t) , intent(in) , value :: nucl_num
|
||||
integer (c_int64_t) , intent(in) , value :: cord_num
|
||||
real (c_double ) , intent(in) , value :: rescale_factor_kappa_en
|
||||
real (c_double ) , intent(in) :: en_distance(nucl_num,elec_num,walk_num)
|
||||
real (c_double ) , intent(out) :: een_rescaled_n(0:cord_num,nucl_num,elec_num,walk_num)
|
||||
|
||||
integer(c_int32_t), external :: qmckl_compute_een_rescaled_n_f
|
||||
info = qmckl_compute_een_rescaled_n_f &
|
||||
(context, &
|
||||
walk_num, &
|
||||
elec_num, &
|
||||
nucl_num, &
|
||||
cord_num, &
|
||||
rescale_factor_kappa_en, &
|
||||
en_distance, &
|
||||
een_rescaled_n)
|
||||
|
||||
end function qmckl_compute_een_rescaled_n
|
||||
#+end_src
|
||||
|
||||
*** Test
|
||||
|
||||
#+begin_src python :results output :exports none :noweb yes
|
||||
import numpy as np
|
||||
|
||||
<<jastrow_data>>
|
||||
|
||||
elec_coord = np.array(elec_coord)[0]
|
||||
nucl_coord = np.array(nucl_coord)
|
||||
elnuc_dist = np.zeros(shape=(elec_num, nucl_num),dtype=float)
|
||||
for i in range(elec_num):
|
||||
for a in range(nucl_num):
|
||||
elnuc_dist[i, a] = np.linalg.norm(elec_coord[i] - nucl_coord[:,a])
|
||||
|
||||
kappa = 1.0
|
||||
|
||||
een_rescaled_n = np.zeros(shape=(nucl_num, elec_num, cord_num + 1), dtype=float)
|
||||
een_rescaled_n[:,:,0] = 1.0
|
||||
|
||||
for a in range(nucl_num):
|
||||
for i in range(elec_num):
|
||||
een_rescaled_n[a, i, 1] = np.exp(-kappa * elnuc_dist[i, a])
|
||||
|
||||
for l in range(2,cord_num+1):
|
||||
for a in range(nucl_num):
|
||||
for i in range(elec_num):
|
||||
een_rescaled_n[a, i, l] = een_rescaled_n[a, i, l - 1] * een_rescaled_n[a, i, 1]
|
||||
|
||||
print(" een_rescaled_n[0, 2, 1] = ",een_rescaled_n[0, 2, 1])
|
||||
print(" een_rescaled_n[0, 3, 1] = ",een_rescaled_n[0, 3, 1])
|
||||
print(" een_rescaled_n[0, 4, 1] = ",een_rescaled_n[0, 4, 1])
|
||||
print(" een_rescaled_n[1, 3, 2] = ",een_rescaled_n[1, 3, 2])
|
||||
print(" een_rescaled_n[1, 4, 2] = ",een_rescaled_n[1, 4, 2])
|
||||
print(" een_rescaled_n[1, 5, 2] = ",een_rescaled_n[1, 5, 2])
|
||||
#+end_src
|
||||
|
||||
#+RESULTS:
|
||||
: een_rescaled_n[0, 2, 1] = 0.10612983920006765
|
||||
: een_rescaled_n[0, 3, 1] = 0.135652809635553
|
||||
: een_rescaled_n[0, 4, 1] = 0.023391817607642338
|
||||
: een_rescaled_n[1, 3, 2] = 0.880957224822116
|
||||
: een_rescaled_n[1, 4, 2] = 0.027185942659395074
|
||||
: een_rescaled_n[1, 5, 2] = 0.01343938025140174
|
||||
|
||||
#+begin_src c :tangle (eval c_test)
|
||||
assert(qmckl_electron_provided(context));
|
||||
|
||||
double een_rescaled_n[walk_num][elec_num][nucl_num][(cord_num + 1)];
|
||||
rc = qmckl_get_jastrow_een_rescaled_n(context, &(een_rescaled_n[0][0][0][0]));
|
||||
|
||||
// value of (0,2,1)
|
||||
assert(fabs(een_rescaled_n[0][2][0][1]-0.10612983920006765) < 1.e-12);
|
||||
assert(fabs(een_rescaled_n[0][3][0][1]-0.135652809635553) < 1.e-12);
|
||||
assert(fabs(een_rescaled_n[0][4][0][1]-0.023391817607642338) < 1.e-12);
|
||||
assert(fabs(een_rescaled_n[0][3][1][2]-0.880957224822116) < 1.e-12);
|
||||
assert(fabs(een_rescaled_n[0][4][1][2]-0.027185942659395074) < 1.e-12);
|
||||
assert(fabs(een_rescaled_n[0][5][1][2]-0.01343938025140174) < 1.e-12);
|
||||
|
||||
#+end_src
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user