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

Factor ee implemented. #22

Browse Source
This commit is contained in:
vijay gopal chilkuri 2021-07-05 18:32:05 +05:30
parent ad824d0f71
commit 360587ef36
1 changed files with 387 additions and 21 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

406
org/qmckl_jastrow.org
View File

@ -81,26 +81,26 @@ int main() {
| ~double~ | ~aord_vector[aord_num + 1][type_nucl_num]~ | in | Order of a polynomial coefficients | | ~double~ | ~aord_vector[aord_num + 1][type_nucl_num]~ | in | Order of a polynomial coefficients |
| ~double~ | ~bord_vector[bord_num + 1]~ | in | Order of b polynomial coefficients | | ~double~ | ~bord_vector[bord_num + 1]~ | in | Order of b polynomial coefficients |
| ~double~ | ~cord_vector[cord_num][type_nucl_num]~ | in | Order of c polynomial coefficients | | ~double~ | ~cord_vector[cord_num][type_nucl_num]~ | in | Order of c polynomial coefficients |
| ~double~ | ~factor_ee~ | out | Jastrow factor: electron-electron part | | ~double~ | ~factor_ee[walk_num]~ | out | Jastrow factor: electron-electron part |
| ~double~ | ~factor_ee_date~ | out | Jastrow factor: electron-electron part | | ~double~ | ~factor_ee_date~ | out | Jastrow factor: electron-electron part |
| ~double~ | ~factor_en~ | out | Jastrow factor: electron-nucleus part | | ~double~ | ~factor_en[walk_num]~ | out | Jastrow factor: electron-nucleus part |
| ~double~ | ~factor_en_date~ | out | Jastrow factor: electron-nucleus part | | ~double~ | ~factor_en_date~ | out | Jastrow factor: electron-nucleus part |
| ~double~ | ~factor_een~ | out | Jastrow factor: electron-electron-nucleus part | | ~double~ | ~factor_een[walk_num]~ | out | Jastrow factor: electron-electron-nucleus part |
| ~double~ | ~factor_een_date~ | out | Jastrow factor: electron-electron-nucleus part | | ~double~ | ~factor_een_date~ | out | Jastrow factor: electron-electron-nucleus part |
| ~double~ | ~factor_ee_deriv_e[4][nelec]~ | out | Derivative of the Jastrow factor: electron-electron-nucleus part | | ~double~ | ~factor_ee_deriv_e[4][nelec][walk_num]~ | out | Derivative of the Jastrow factor: electron-electron-nucleus part |
| ~double~ | ~factor_en_deriv_e[4][nelec]~ | out | Derivative of the Jastrow factor: electron-electron-nucleus part | | ~double~ | ~factor_en_deriv_e[4][nelec][walk_num]~ | out | Derivative of the Jastrow factor: electron-electron-nucleus part |
| ~double~ | ~factor_een_deriv_e[4][nelec]~ | out | Derivative of the Jastrow factor: electron-electron-nucleus part | | ~double~ | ~factor_een_deriv_e[4][nelec][walk_num]~ | out | Derivative of the Jastrow factor: electron-electron-nucleus part |
computed data: computed data:
|-----------+---------------------------------------------------+-------------------------------------------------| |-----------+-------------------------------------------------------------+-------------------------------------------------|
| ~int64_t~ | ~dim_cord_vec~ | Number of unique C coefficients | | ~int64_t~ | ~dim_cord_vec~ | Number of unique C coefficients |
| ~double~ | ~asymp_jasb[2]~ | Asymptotic component | | ~double~ | ~asymp_jasb[2]~ | Asymptotic component |
| ~int64_t~ | ~asymp_jasb_date~ | Asymptotic component | | ~int64_t~ | ~asymp_jasb_date~ | Asymptotic component |
| ~double~ | ~coord_vect_full[dim_cord_vec][nucl_num]~ | vector of non-zero coefficients | | ~double~ | ~coord_vect_full[dim_cord_vec][nucl_num]~ | vector of non-zero coefficients |
| ~int64_t~ | ~lkpm_of_cindex[4][dim_cord_vec]~ | Transform l,k,p, and m into consecutive indices | | ~int64_t~ | ~lkpm_of_cindex[4][dim_cord_vec]~ | Transform l,k,p, and m into consecutive indices |
| ~double~ | ~tmp_c[elec_num][nucl_num][ncord + 1][ncord]~ | vector of non-zero coefficients | | ~double~ | ~tmp_c[elec_num][nucl_num][ncord + 1][ncord][walk_num]~ | vector of non-zero coefficients |
| ~double~ | ~dtmp_c[elec_num][4][nucl_num][ncord + 1][ncord]~ | vector of non-zero coefficients | | ~double~ | ~dtmp_c[elec_num][4][nucl_num][ncord + 1][ncord][walk_num]~ | vector of non-zero coefficients |
For H2O we have the following data: For H2O we have the following data:
@ -108,6 +108,68 @@ int main() {
#+BEGIN_SRC python :results output #+BEGIN_SRC python :results output
import numpy as np import numpy as np
elec_num = 10
nucl_num = 2
up_num = 5
down_num = 5
nucl_coord = [ [0.000000, 0.000000 ],
[0.000000, 0.000000 ],
[2.059801, 2.059801 ] ]
elec_coord = [[[-0.250655104764153 , 0.503070975550133 , -0.166554344502303],
[-0.587812193472177 , -0.128751981129274 , 0.187773606533075],
[ 1.61335569047166 , -0.615556732874863 , -1.43165470979934 ],
[-4.901239896295210E-003 , -1.120440036458986E-002 , 1.99761909330422 ],
[ 0.766647499681200 , -0.293515395797937 , 3.66454589201239 ],
[-0.127732483187947 , -0.138975497694196 , -8.669850480215846E-002],
[-0.232271834949124 , -1.059321673434182E-002 , -0.504862241464867],
[ 1.09360863531826 , -2.036103063808752E-003 , -2.702796910818986E-002],
[-0.108090166832043 , 0.189161729653261 , 2.15398313919894],
[ 0.397978144318712 , -0.254277292595981 , 2.54553335476344]]];
ee_distance_rescaled = [
[ 0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
0.000000000000000E+000],
[ 0.550227800352402 ,0.000000000000000E+000 ,0.000000000000000E+000,
0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
0.000000000000000E+000],
[ 0.919155060185168 ,0.937695909123175 ,0.000000000000000E+000,
0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
0.000000000000000E+000],
[ 0.893325429242815 ,0.851181978173561 ,0.978501685226877 ,
0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
0.000000000000000E+000],
[ 0.982457268305353 ,0.976125002619471 ,0.994349933143149 ,
0.844077311588328 ,0.000000000000000E+000 ,0.000000000000000E+000,
0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
0.000000000000000E+000],
[ 0.482407528408731 ,0.414816073699124 ,0.894716035479343 ,
0.876540187084407 ,0.978921170036895 ,0.000000000000000E+000,
0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
0.000000000000000E+000],
[ 0.459541909660400 ,0.545007215761510 ,0.883752955884551 ,
0.918958134888791 ,0.986386936267237 ,0.362209822236419 ,
0.000000000000000E+000 ,0.000000000000000E+000 ,0.000000000000000E+000,
0.000000000000000E+000],
[ 0.763732576854455 ,0.817282762358449 ,0.801802919535959 ,
0.900089095449775 ,0.975704636491453 ,0.707836537586060 ,
0.755705808346586 ,0.000000000000000E+000 ,0.000000000000000E+000,
0.000000000000000E+000],
[ 0.904249454052971 ,0.871097965261373 ,0.982717262706270 ,
0.239901207363622 ,0.836519456769083 ,0.896135326270534 ,
0.930694340243023 ,0.917708540815567 ,0.000000000000000E+000,
0.000000000000000E+000],
[ 0.944400908070716 ,0.922589018494961 ,0.984615718580670 ,
0.514328661540623 ,0.692362267147064 ,0.931894098453677 ,
0.956034127544344 ,0.931221472309472 ,0.540903688625053 ,
0.000000000000000E+000]]
type_nucl_num = 1 type_nucl_num = 1
aord_num = 5 aord_num = 5
bord_num = 5 bord_num = 5
@ -167,6 +229,8 @@ lkpm_of_cindex = [[1 , 1 , 2 , 0],
[3 , 0 , 5 , 1], [3 , 0 , 5 , 1],
[1 , 0 , 5 , 2]] [1 , 0 , 5 , 2]]
kappa = 1.0
kappa_inv = 1.0/kappa
#+END_SRC #+END_SRC
#+RESULTS: jastrow_data #+RESULTS: jastrow_data
@ -976,7 +1040,7 @@ assert(qmckl_nucleus_provided(context));
compute. If it is the case, then the data is recomputed and the compute. If it is the case, then the data is recomputed and the
current date is stored. current date is stored.
** Asymptotic component for \(J_{ee}\) ** Asymptotic component for \(f_{ee}\)
Calculate the asymptotic component ~asymp_jasb~ to be substracted from the final Calculate the asymptotic component ~asymp_jasb~ to be substracted from the final
electron-electron jastrow factor \(f_{ee}\). The asymptotic componenet is calculated electron-electron jastrow factor \(f_{ee}\). The asymptotic componenet is calculated
@ -1107,13 +1171,11 @@ integer function qmckl_compute_asymp_jasb_f(context, bord_num, bord_vector, resc
info = QMCKL_INVALID_CONTEXT info = QMCKL_INVALID_CONTEXT
return return
endif endif
print *,"In Compute 1", asym_one
if (bord_num <= 0) then if (bord_num <= 0) then
info = QMCKL_INVALID_ARG_2 info = QMCKL_INVALID_ARG_2
return return
endif endif
print *,"In Compute 2", asym_one
asym_one = bord_vector(1) * kappa_inv / (1.0d0 + bord_vector(2) * kappa_inv) asym_one = bord_vector(1) * kappa_inv / (1.0d0 + bord_vector(2) * kappa_inv)
asymp_jasb(:) = (/asym_one, 0.5d0 * asym_one/) asymp_jasb(:) = (/asym_one, 0.5d0 * asym_one/)
@ -1167,14 +1229,12 @@ end function qmckl_compute_asymp_jasb_f
#+end_src #+end_src
*** Test *** Test
#+name: asymp_jasb
#+begin_src python :results output :exports none :noweb yes #+begin_src python :results output :exports none :noweb yes
import numpy as np import numpy as np
<<jastrow_data>> <<jastrow_data>>
kappa = 1.0
kappa_inv = 1.0
asym_one = bord_vector[0] * kappa_inv / (1.0 + bord_vector[1]*kappa_inv) asym_one = bord_vector[0] * kappa_inv / (1.0 + bord_vector[1]*kappa_inv)
asymp_jasb = np.array([asym_one, 0.5 * asym_one]) asymp_jasb = np.array([asym_one, 0.5 * asym_one])
@ -1189,6 +1249,11 @@ print("asymp_jasb[0] : ", asymp_jasb[0])
print("asymp_jasb[1] : ", asymp_jasb[1]) print("asymp_jasb[1] : ", asymp_jasb[1])
#+end_src #+end_src
#+RESULTS: asymp_jasb
: asym_one : 0.6634291325000664
: asymp_jasb[0] : 1.043287918508297
: asymp_jasb[1] : 0.7115733522582638
#+RESULTS: #+RESULTS:
: asym_one : 0.43340325572525706 : asym_one : 0.43340325572525706
: asymp_jasb[0] : 0.5323750557252571 : asymp_jasb[0] : 0.5323750557252571
@ -1232,6 +1297,307 @@ assert(fabs(asymp_jasb[1]-0.31567342786262853) < 1.e-12);
#+end_src #+end_src
** Electron-electron component \(f_{ee}\)
Calculate the asymptotic component ~asymp_jasb~ to be substracted from the final
electron-electron jastrow factor \(f_{ee}\). The asymptotic componenet is calculated
via the ~bord_vector~ and the electron-electron rescale factor ~rescale_factor_kappa~.
\[
f_{ee} = \sum_{i,j<i} \left\{ \frac{ \eta B_0 C_{ij}}{1 - B_1 C_{ij}} - J_{asymp} + \sum^{nord}_{k}B_k C_{ij}^k \right\}
\]
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_jastrow_factor_ee(qmckl_context context, double* const factor_ee);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_jastrow_factor_ee(qmckl_context context, double* const factor_ee)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_exit_code rc;
rc = qmckl_provide_factor_ee(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL);
memcpy(factor_ee, ctx->jastrow.factor_ee, ctx->electron.walk_num*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_factor_ee(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_provide_factor_ee(qmckl_context context)
{
qmckl_exit_code rc;
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 rescaled distance is provided */
rc = qmckl_provide_ee_distance_rescaled(context);
if(rc != QMCKL_SUCCESS) return rc;
/* Compute if necessary */
if (ctx->date > ctx->jastrow.factor_ee_date) {
/* Allocate array */
if (ctx->jastrow.factor_ee == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->electron.walk_num * sizeof(double);
double* factor_ee = (double*) qmckl_malloc(context, mem_info);
if (factor_ee == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_provide_factor_ee",
NULL);
}
ctx->jastrow.factor_ee = factor_ee;
}
qmckl_exit_code rc =
qmckl_compute_factor_ee(context,
ctx->electron.walk_num,
ctx->electron.num,
ctx->electron.up_num,
ctx->jastrow.bord_num,
ctx->jastrow.bord_vector,
ctx->electron.ee_distance_rescaled,
ctx->jastrow.asymp_jasb,
ctx->jastrow.factor_ee);
if (rc != QMCKL_SUCCESS) {
return rc;
}
ctx->jastrow.factor_ee_date = ctx->date;
}
return QMCKL_SUCCESS;
}
#+end_src
*** Compute
:PROPERTIES:
:Name: qmckl_compute_factor_ee
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
#+NAME: qmckl_factor_ee_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 | up_num | in | Number of alpha electrons |
| int64_t | bord_num | in | Number of coefficients |
| double | bord_vector[bord_num + 1] | in | List of coefficients |
| double | ee_distance_rescaled[walk_num][elec_num][elec_num] | in | Electron-electron distances |
| double | asymp_jasb[2] | in | Electron-electron distances |
| double | factor_ee[walk_num] | out | Electron-electron distances |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_factor_ee_f(context, walk_num, elec_num, up_num, bord_num, &
bord_vector, ee_distance_rescaled, asymp_jasb, factor_ee) &
result(info)
use qmckl
implicit none
integer(qmckl_context), intent(in) :: context
integer*8 , intent(in) :: walk_num, elec_num, bord_num, up_num
double precision , intent(in) :: bord_vector(bord_num)
double precision , intent(in) :: ee_distance_rescaled(walk_num, elec_num, elec_num)
double precision , intent(in) :: asymp_jasb(2)
double precision , intent(out) :: factor_ee(walk_num)
integer*8 :: i, j, p, ipar, nw
double precision :: pow_ser, x, spin_fact, power_ser
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 (bord_num <= 0) then
info = QMCKL_INVALID_ARG_4
return
endif
factor_ee = 0.0d0
do nw =1, walk_num
do j = 1, elec_num
do i = 1, j - 1
x = ee_distance_rescaled(nw,i,j)
power_ser = 0.0d0
spin_fact = 1.0d0
ipar = 1
do p = 2, bord_num
x = x * ee_distance_rescaled(nw,i,j)
power_ser = power_ser + bord_vector(p + 1) * x
end do
if(j .LE. up_num .OR. i .GT. up_num) then
spin_fact = 0.5d0
ipar = 2
endif
factor_ee(nw) = factor_ee(nw) + spin_fact * bord_vector(1) * &
ee_distance_rescaled(nw,i,j) / &
(1.0d0 + bord_vector(2) * &
ee_distance_rescaled(nw,i,j)) &
-asymp_jasb(ipar) + power_ser
end do
end do
end do
end function qmckl_compute_factor_ee_f
#+end_src
#+CALL: generate_c_header(table=qmckl_factor_ee_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
qmckl_exit_code qmckl_compute_factor_ee (
const qmckl_context context,
const int64_t walk_num,
const int64_t elec_num,
const int64_t up_num,
const int64_t bord_num,
const double* bord_vector,
const double* ee_distance_rescaled,
const double* asymp_jasb,
double* const factor_ee );
#+end_src
#+CALL: generate_c_interface(table=qmckl_factor_ee_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_factor_ee &
(context, &
walk_num, &
elec_num, &
up_num, &
bord_num, &
bord_vector, &
ee_distance_rescaled, &
asymp_jasb, &
factor_ee) &
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 :: up_num
integer (c_int64_t) , intent(in) , value :: bord_num
real (c_double ) , intent(in) :: bord_vector(bord_num + 1)
real (c_double ) , intent(in) :: ee_distance_rescaled(elec_num,elec_num,walk_num)
real (c_double ) , intent(in) :: asymp_jasb(2)
real (c_double ) , intent(out) :: factor_ee(walk_num)
integer(c_int32_t), external :: qmckl_compute_factor_ee_f
info = qmckl_compute_factor_ee_f &
(context, &
walk_num, &
elec_num, &
up_num, &
bord_num, &
bord_vector, &
ee_distance_rescaled, &
asymp_jasb, &
factor_ee)
end function qmckl_compute_factor_ee
#+end_src
*** Test
#+begin_src python :results output :exports none :noweb yes
import numpy as np
<<jastrow_data>>
<<asymp_jasb>>
factor_ee = 0.0
for i in range(0,elec_num):
for j in range(0,i):
x = ee_distance_rescaled[i][j]
pow_ser = 0.0
spin_fact = 1.0
ipar = 0
for p in range(1,bord_num):
x = x * ee_distance_rescaled[i][j]
pow_ser = pow_ser + bord_vector[p + 1] * x
if(i < up_num or j >= up_num):
spin_fact = 0.5
ipar = 1
factor_ee = factor_ee + spin_fact * bord_vector[0] * ee_distance_rescaled[i][j] \
/ (1.0 + bord_vector[1] * ee_distance_rescaled[i][j]) \
- asymp_jasb[ipar] + pow_ser
print("factor_ee :",factor_ee)
#+end_src
#+RESULTS:
: asym_one : 0.43340325572525706
: asymp_jasb[0] : 0.5323750557252571
: asymp_jasb[1] : 0.31567342786262853
: factor_ee : -4.282760865958113
#+begin_src c :tangle (eval c_test)
/* Check if Jastrow is properly initialized */
assert(qmckl_jastrow_provided(context));
double factor_ee[walk_num];
rc = qmckl_get_jastrow_factor_ee(context, factor_ee);
// calculate factor_ee
assert(fabs(factor_ee[0]+4.282760865958113) < 1.e-12);
#+end_src
* End of files :noexport: * End of files :noexport:
#+begin_src c :tangle (eval h_private_type) #+begin_src c :tangle (eval h_private_type)