TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-12-22 20:36:01 +01:00
Code Issues Projects Releases Wiki Activity
24e3f8dd11
qmckl/org/qmckl_local_energy.org

1739 lines
58 KiB
Org Mode
Raw Blame History

#+TITLE: Local Energy
#+SETUPFILE: ../tools/theme.setup
#+INCLUDE: ../tools/lib.org
Here we calculate the final expectation value of the
local energy \[E_L\] as the sum of the kinetic energy and potential energy.
\[
E_L = KE + PE
\]
Where the kinetic energy is given as:
\[
KE = -\frac{1}{2}\frac{\bigtriangleup \Psi}{\Psi}
\]
The laplacian of the wavefunction in the single-determinant
case is given as follows:
\[
\frac{\bigtriangleup \Psi(r)}{\Psi(r)} = \sum_{j=1}^{N_e} \bigtriangleup \Phi_j(r_i) D_{ji}^{-1}(r)
\]
The potential energy is the sum of all the following terms
\[
PE = \mathcal{V}_{ee} + \mathcal{V}_{en} + \mathcal{V}_{nn}
\]
The potential for is calculated as the sum of single electron
contributions.
\[
\mathcal{V}_{ee} = \sum_{i=1}^{N_e}\sum_{j<i} \frac{1}{r_{ij}}
\]
\[
\mathcal{V}_{en} = - \sum_{i=1}^{N_e}\sum_{A=1}^{N_n}\frac{Z_A}{r_{iA}}
\]
\[
\mathcal{V}_{nn} = \sum_{A=1}^{N_n}\sum_{B<A}\frac{Z_A Z_B}{r_{AB}}
\]
The remaining quantities that are required for the calculation of
a single Monte-Carlo step are as follows:
1. Drift Vector - \[F(x)\]
\[
F(x) = 2\frac{\nabla \Psi(r)}{\Psi(r)}
\]
2. Diffusion move - \[y\]
\[
y = x + D F(x) \delta t + \chi
\]
Where \[\chi\] is a random number with gaussian distribution centered at 0
and width of \[2D\delta t\]. Here \[D\] is the drift parameter.
3. Acceptance probability - \[min\left[1, q(y,x)\right]\]
\[
q(y,x) = Exp\left[\frac{\delta r}{2} ( F(x) + F(y)) - \frac{D\delta t}{4}\left( F(x)^2 - F(y)^2\right)\right]
\]
With these quantities, a single determinant VMC simulation can be carried out.
* Headers :noexport:
#+begin_src elisp :noexport :results none
(org-babel-lob-ingest "../tools/lib.org")
#+end_src
#+begin_src c :tangle (eval h_private_type)
#ifndef QMCKL_LOCAL_ENERGY_HPT
#define QMCKL_LOCAL_ENERGY_HPT
#include <stdbool.h>
#+end_src
#+begin_src c :tangle (eval c_test) :noweb yes
#include "qmckl.h"
#include "assert.h"
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <stdio.h>
#include <math.h>
#include "chbrclf.h"
#include "qmckl_ao_private_func.h"
#include "qmckl_mo_private_func.h"
#include "qmckl_determinant_private_func.h"
int main() {
qmckl_context context;
context = qmckl_context_create();
qmckl_exit_code rc;
#+end_src
#+begin_src c :tangle (eval c)
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#ifdef HAVE_STDINT_H
#include <stdint.h>
#elif HAVE_INTTYPES_H
#include <inttypes.h>
#endif
#include <stdlib.h>
#include <string.h>
#include <stdbool.h>
#include <assert.h>
#include "qmckl.h"
#include "qmckl_context_private_type.h"
#include "qmckl_memory_private_type.h"
#include "qmckl_memory_private_func.h"
#include "qmckl_ao_private_type.h"
#include "qmckl_ao_private_func.h"
#include "qmckl_mo_private_type.h"
#include "qmckl_mo_private_func.h"
#include "qmckl_determinant_private_type.h"
#include "qmckl_determinant_private_func.h"
#+end_src
* Context
The following arrays are stored in the context:
| |
Computed data:
|--------------+---------------------------+----------------------------|
| ~e_kin~ | ~[walk_num]~ | total kinetic energy |
| ~e_pot~ | ~[walk_num]~ | total potential energy |
| ~e_local~ | ~[walk_num]~ | local energy |
| ~r_drift~ | ~[3][walk_num][elec_num]~ | The drift vector |
| ~y_move~ | ~[3][walk_num]~ | The diffusion move |
| ~accep_prob~ | ~[walk_num]~ | The acceptance probability |
|--------------+---------------------------+----------------------------|
** Data structure
#+begin_src c :comments org :tangle (eval h_private_type)
typedef struct qmckl_local_energy_struct {
double * e_kin;
double * e_pot;
double * e_local;
double * accep_prob;
double * r_drift;
double * y_move;
uint64_t e_kin_date;
uint64_t e_pot_date;
uint64_t e_local_date;
uint64_t accep_prob_date;
uint64_t r_drift_date;
uint64_t y_move_date;
int32_t uninitialized;
bool provided;
} qmckl_local_energy_struct;
#+end_src
The ~uninitialized~ integer contains one bit set to one for each
initialization function which has not been called. It becomes equal
to zero after all initialization functions have been called. The
struct is then initialized and ~provided == true~.
Some values are initialized by default, and are not concerned by
this mechanism.
** Access functions
When all the data for the local energy have been provided, the following
function returns ~true~.
#+begin_src c :comments org :tangle (eval h_func)
bool qmckl_local_energy_provided (const qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
bool qmckl_local_energy_provided(const qmckl_context context) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return false;
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
qmckl_exit_code rc;
if(!qmckl_electron_provided(context)) return false;
if(!qmckl_nucleus_provided(context)) return false;
rc = qmckl_provide_ao_basis_ao_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_mo_basis_mo_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
ctx->local_energy.provided = (ctx->local_energy.uninitialized == 0);
return ctx->local_energy.provided;
}
#+end_src
* Computation
** Kinetic energy
:PROPERTIES:
:Name: qmckl_compute_kinetic_energy
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
Where the kinetic energy is given as:
\[
KE = -\frac{1}{2}\frac{\bigtriangleup \Psi}{\Psi}
\]
The laplacian of the wavefunction in the single-determinant
case is given as follows:
\[
\frac{\bigtriangleup \Psi(r)}{\Psi(r)} = \sum_{j=1}^{N_e} \bigtriangleup \Phi_j(r_i) D_{ji}^{-1}(r)
\]
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_kinetic_energy(qmckl_context context, double* const kinetic_energy);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_kinetic_energy(qmckl_context context, double * const kinetic_energy) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_exit_code rc;
if(!qmckl_electron_provided(context)) return QMCKL_NOT_PROVIDED;
if(!qmckl_nucleus_provided(context)) return QMCKL_NOT_PROVIDED;
rc = qmckl_provide_ao_basis_ao_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_mo_basis_mo_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_kinetic_energy(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->electron.walker.num * sizeof(double);
memcpy(kinetic_energy, ctx->local_energy.e_kin, sze);
return QMCKL_SUCCESS;
}
#+end_src
*** Provide
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_kinetic_energy(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_provide_kinetic_energy(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*) context;
assert (ctx != NULL);
if(!(ctx->nucleus.provided)) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_electron",
NULL);
}
if(!(ctx->electron.provided)) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_electron",
NULL);
}
if (!ctx->ao_basis.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_ao_basis",
NULL);
}
if (!ctx->mo_basis.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_mo_basis",
NULL);
}
if (!ctx->det.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_mo_basis",
NULL);
}
rc = qmckl_provide_det_inv_matrix_alpha(context);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_det_inv_matrix_alpha",
NULL);
}
rc = qmckl_provide_det_inv_matrix_beta(context);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_det_inv_matrix_beta",
NULL);
}
/* Compute if necessary */
if (ctx->electron.walker.point.date > ctx->local_energy.e_kin_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->local_energy.e_kin);
ctx->local_energy.e_kin = NULL;
}
/* Allocate array */
if (ctx->local_energy.e_kin == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->electron.walker.num * sizeof(double);
double* e_kin = (double*) qmckl_malloc(context, mem_info);
if (e_kin == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_e_kin",
NULL);
}
ctx->local_energy.e_kin = e_kin;
}
if (ctx->det.type == 'G') {
rc = qmckl_compute_kinetic_energy(context,
ctx->electron.walker.num,
ctx->det.det_num_alpha,
ctx->det.det_num_beta,
ctx->electron.up_num,
ctx->electron.down_num,
ctx->electron.num,
ctx->det.mo_index_alpha,
ctx->det.mo_index_beta,
ctx->mo_basis.mo_num,
ctx->mo_basis.mo_vgl,
ctx->det.det_value_alpha,
ctx->det.det_value_beta,
ctx->det.det_inv_matrix_alpha,
ctx->det.det_inv_matrix_beta,
ctx->local_energy.e_kin);
} else {
return qmckl_failwith( context,
QMCKL_FAILURE,
"compute_kinetic_energy",
"Not yet implemented");
}
if (rc != QMCKL_SUCCESS) {
return rc;
}
ctx->local_energy.e_kin_date = ctx->date;
}
return QMCKL_SUCCESS;
}
#+end_src
*** Compute kinetic enregy
:PROPERTIES:
:Name: qmckl_compute_kinetic_energy
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
#+NAME: qmckl_compute_kinetic_energy_args
| ~qmckl_context~ | ~context~ | in | Global state |
| ~int64_t~ | ~walk_num~ | in | Number of walkers |
| ~int64_t~ | ~det_num_alpha~ | in | Number of determinants |
| ~int64_t~ | ~det_num_beta~ | in | Number of determinants |
| ~int64_t~ | ~alpha_num~ | in | Number of electrons |
| ~int64_t~ | ~beta_num~ | in | Number of electrons |
| ~int64_t~ | ~elec_num~ | in | Number of electrons |
| ~int64_t~ | ~mo_index_alpha[det_num_alpha][walk_num][alpha_num]~ | in | MO indices for electrons |
| ~int64_t~ | ~mo_index_beta[det_num_beta][walk_num][beta_num]~ | in | MO indices for electrons |
| ~int64_t~ | ~mo_num~ | in | Number of MOs |
| ~double~ | ~mo_vgl[5][elec_num][mo_num]~ | in | Value, gradients and Laplacian of the MOs |
| ~double~ | ~det_value_alpha[det_num_alpha][walk_num]~ | in | Det of wavefunction |
| ~double~ | ~det_value_beta[det_num_beta][walk_num]~ | in | Det of wavefunction |
| ~double~ | ~det_inv_matrix_alpha[det_num_alpha][walk_num][alpha_num][alpha_num]~ | in | Value, gradients and Laplacian of the Det |
| ~double~ | ~det_inv_matrix_beta[det_num_beta][walk_num][beta_num][beta_num]~ | in | Value, gradients and Laplacian of the Det |
| ~double~ | ~e_kin[walk_num]~ | out | Kinetic energy |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_kinetic_energy_f(context, walk_num, &
det_num_alpha, det_num_beta, alpha_num, beta_num, elec_num, mo_index_alpha, mo_index_beta, &
mo_num, mo_vgl, det_value_alpha, det_value_beta, det_inv_matrix_alpha, det_inv_matrix_beta, e_kin) &
result(info)
use qmckl
implicit none
integer(qmckl_context) , intent(in) :: context
integer*8, intent(in) :: walk_num
integer*8, intent(in) :: det_num_alpha
integer*8, intent(in) :: det_num_beta
integer*8, intent(in) :: alpha_num
integer*8, intent(in) :: beta_num
integer*8, intent(in) :: elec_num
integer*8, intent(in) :: mo_num
integer*8, intent(in) :: mo_index_alpha(alpha_num, walk_num, det_num_alpha)
integer*8, intent(in) :: mo_index_beta(beta_num, walk_num, det_num_beta)
double precision, intent(in) :: mo_vgl(mo_num, elec_num, 5)
double precision, intent(in) :: det_value_alpha(walk_num, det_num_alpha)
double precision, intent(in) :: det_value_beta(walk_num, det_num_beta)
double precision, intent(in) :: det_inv_matrix_alpha(alpha_num, alpha_num, walk_num, det_num_alpha)
double precision, intent(in) :: det_inv_matrix_beta(beta_num, beta_num, walk_num, det_num_beta)
double precision, intent(inout) :: e_kin(walk_num)
double precision :: tmp_e
integer*8 :: idet, iwalk, ielec, mo_id, imo
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 (alpha_num <= 0) then
info = QMCKL_INVALID_ARG_3
return
endif
if (beta_num < 0) then
info = QMCKL_INVALID_ARG_4
return
endif
if (elec_num <= 0) then
info = QMCKL_INVALID_ARG_5
return
endif
e_kin = 0.0d0
do idet = 1, det_num_alpha
do iwalk = 1, walk_num
! Alpha part
do imo = 1, alpha_num
do ielec = 1, alpha_num
mo_id = mo_index_alpha(imo, iwalk, idet)
e_kin(iwalk) = e_kin(iwalk) - 0.5d0 * det_inv_matrix_alpha(imo, ielec, iwalk, idet) * &
mo_vgl(mo_id, ielec, 5)
end do
end do
! Beta part
do imo = 1, beta_num
do ielec = 1, beta_num
mo_id = mo_index_beta(imo, iwalk, idet)
e_kin(iwalk) = e_kin(iwalk) - 0.5d0 * det_inv_matrix_beta(imo, ielec, iwalk, idet) * &
mo_vgl(mo_id, alpha_num + ielec, 5)
end do
end do
end do
end do
end function qmckl_compute_kinetic_energy_f
#+end_src
#+CALL: generate_c_header(table=qmckl_compute_kinetic_energy_args,rettyp=get_value("CRetType"),fname="qmckl_compute_kinetic_energy"))
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
qmckl_exit_code qmckl_compute_kinetic_energy (
const qmckl_context context,
const int64_t walk_num,
const int64_t det_num_alpha,
const int64_t det_num_beta,
const int64_t alpha_num,
const int64_t beta_num,
const int64_t elec_num,
const int64_t* mo_index_alpha,
const int64_t* mo_index_beta,
const int64_t mo_num,
const double* mo_vgl,
const double* det_value_alpha,
const double* det_value_beta,
const double* det_inv_matrix_alpha,
const double* det_inv_matrix_beta,
double* const e_kin );
#+end_src
#+CALL: generate_c_interface(table=qmckl_compute_kinetic_energy_args,rettyp=get_value("CRetType"),fname="qmckl_compute_kinetic_energy"))
#+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none
integer(c_int32_t) function qmckl_compute_kinetic_energy &
(context, &
walk_num, &
det_num_alpha, &
det_num_beta, &
alpha_num, &
beta_num, &
elec_num, &
mo_index_alpha, &
mo_index_beta, &
mo_num, &
mo_vgl, &
det_value_alpha, &
det_value_beta, &
det_inv_matrix_alpha, &
det_inv_matrix_beta, &
e_kin) &
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 :: det_num_alpha
integer (c_int64_t) , intent(in) , value :: det_num_beta
integer (c_int64_t) , intent(in) , value :: alpha_num
integer (c_int64_t) , intent(in) , value :: beta_num
integer (c_int64_t) , intent(in) , value :: elec_num
integer (c_int64_t) , intent(in) :: mo_index_alpha(alpha_num,walk_num,det_num_alpha)
integer (c_int64_t) , intent(in) :: mo_index_beta(beta_num,walk_num,det_num_beta)
integer (c_int64_t) , intent(in) , value :: mo_num
real (c_double ) , intent(in) :: mo_vgl(mo_num,elec_num,5)
real (c_double ) , intent(in) :: det_value_alpha(walk_num,det_num_alpha)
real (c_double ) , intent(in) :: det_value_beta(walk_num,det_num_beta)
real (c_double ) , intent(in) :: det_inv_matrix_alpha(alpha_num,alpha_num,walk_num,det_num_alpha)
real (c_double ) , intent(in) :: det_inv_matrix_beta(beta_num,beta_num,walk_num,det_num_beta)
real (c_double ) , intent(out) :: e_kin(walk_num)
integer(c_int32_t), external :: qmckl_compute_kinetic_energy_f
info = qmckl_compute_kinetic_energy_f &
(context, &
walk_num, &
det_num_alpha, &
det_num_beta, &
alpha_num, &
beta_num, &
elec_num, &
mo_index_alpha, &
mo_index_beta, &
mo_num, &
mo_vgl, &
det_value_alpha, &
det_value_beta, &
det_inv_matrix_alpha, &
det_inv_matrix_beta, &
e_kin)
end function qmckl_compute_kinetic_energy
#+end_src
*** Test
#+begin_src c :tangle (eval c_test) :exports none
double* elec_coord = &(chbrclf_elec_coord[0][0][0]);
const double* nucl_charge = chbrclf_charge;
const double* nucl_coord = &(chbrclf_nucl_coord[0][0]);
rc = qmckl_set_electron_num (context, chbrclf_elec_up_num, chbrclf_elec_dn_num);
assert (rc == QMCKL_SUCCESS);
assert(qmckl_electron_provided(context));
rc = qmckl_set_electron_coord (context, 'N', chbrclf_walk_num, elec_coord, chbrclf_walk_num*chbrclf_elec_num*3);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_num (context, chbrclf_nucl_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), chbrclf_nucl_num*3);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_nucleus_charge(context, nucl_charge, chbrclf_nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(qmckl_nucleus_provided(context));
const int64_t * nucleus_index = &(chbrclf_basis_nucleus_index[0]);
const int64_t * nucleus_shell_num = &(chbrclf_basis_nucleus_shell_num[0]);
const int32_t * shell_ang_mom = &(chbrclf_basis_shell_ang_mom[0]);
const int64_t * shell_prim_num = &(chbrclf_basis_shell_prim_num[0]);
const int64_t * shell_prim_index = &(chbrclf_basis_shell_prim_index[0]);
const double * shell_factor = &(chbrclf_basis_shell_factor[0]);
const double * exponent = &(chbrclf_basis_exponent[0]);
const double * coefficient = &(chbrclf_basis_coefficient[0]);
const double * prim_factor = &(chbrclf_basis_prim_factor[0]);
const double * ao_factor = &(chbrclf_basis_ao_factor[0]);
const char typ = 'G';
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_type (context, typ);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_shell_num (context, chbrclf_shell_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_prim_num (context, chbrclf_prim_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_nucleus_index (context, nucleus_index, chbrclf_nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_nucleus_shell_num (context, nucleus_shell_num, chbrclf_nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_shell_ang_mom (context, shell_ang_mom, chbrclf_shell_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_shell_factor (context, shell_factor, chbrclf_shell_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_shell_prim_num (context, shell_prim_num, chbrclf_shell_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_shell_prim_index (context, shell_prim_index, chbrclf_shell_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_exponent (context, exponent, chbrclf_prim_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_coefficient (context, coefficient, chbrclf_prim_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_ao_basis_provided(context));
rc = qmckl_set_ao_basis_prim_factor (context, prim_factor, chbrclf_prim_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_ao_basis_ao_num(context, chbrclf_ao_num);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_ao_basis_ao_factor (context, ao_factor, chbrclf_ao_num);
assert(rc == QMCKL_SUCCESS);
assert(qmckl_ao_basis_provided(context));
double ao_vgl[chbrclf_walk_num*chbrclf_elec_num][5][chbrclf_ao_num];
rc = qmckl_get_ao_basis_ao_vgl(context, &(ao_vgl[0][0][0]),
(int64_t) 5*chbrclf_walk_num*chbrclf_elec_num*chbrclf_ao_num);
assert (rc == QMCKL_SUCCESS);
/* Set up MO data */
const int64_t mo_num = chbrclf_mo_num;
rc = qmckl_set_mo_basis_mo_num(context, mo_num);
assert (rc == QMCKL_SUCCESS);
const double * mo_coefficient = &(chbrclf_mo_coef[0]);
rc = qmckl_set_mo_basis_coefficient(context, mo_coefficient, chbrclf_ao_num*chbrclf_mo_num);
assert (rc == QMCKL_SUCCESS);
assert(qmckl_mo_basis_provided(context));
double mo_vgl[chbrclf_walk_num*chbrclf_elec_num][5][chbrclf_mo_num];
rc = qmckl_get_mo_basis_mo_vgl(context, &(mo_vgl[0][0][0]), 5*chbrclf_walk_num*chbrclf_elec_num*chbrclf_mo_num);
assert (rc == QMCKL_SUCCESS);
/* Set up determinant data */
#define det_num_alpha 1
#define det_num_beta 1
int64_t mo_index_alpha[det_num_alpha][chbrclf_walk_num][chbrclf_elec_up_num];
int64_t mo_index_beta[det_num_alpha][chbrclf_walk_num][chbrclf_elec_dn_num];
int i, j, k;
for(k = 0; k < det_num_alpha; ++k)
for(i = 0; i < chbrclf_walk_num; ++i)
for(j = 0; j < chbrclf_elec_up_num; ++j)
mo_index_alpha[k][i][j] = j + 1;
for(k = 0; k < det_num_beta; ++k)
for(i = 0; i < chbrclf_walk_num; ++i)
for(j = 0; j < chbrclf_elec_dn_num; ++j)
mo_index_beta[k][i][j] = j + 1;
rc = qmckl_set_determinant_type (context, typ);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_set_determinant_det_num_alpha (context, det_num_alpha);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_set_determinant_det_num_beta (context, det_num_beta);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_set_determinant_mo_index_alpha (context, &(mo_index_alpha[0][0][0]), det_num_alpha*chbrclf_walk_num*chbrclf_elec_up_num);
assert (rc == QMCKL_SUCCESS);
rc = qmckl_set_determinant_mo_index_beta (context, &(mo_index_beta[0][0][0]),det_num_beta*chbrclf_walk_num*chbrclf_elec_dn_num);
assert (rc == QMCKL_SUCCESS);
// Get alpha determinant
double det_vgl_alpha[det_num_alpha][chbrclf_walk_num][5][chbrclf_elec_up_num][chbrclf_elec_up_num];
double det_vgl_beta[det_num_beta][chbrclf_walk_num][5][chbrclf_elec_dn_num][chbrclf_elec_dn_num];
rc = qmckl_get_det_vgl_alpha(context, &(det_vgl_alpha[0][0][0][0][0]));
assert (rc == QMCKL_SUCCESS);
rc = qmckl_get_det_vgl_beta(context, &(det_vgl_beta[0][0][0][0][0]));
assert (rc == QMCKL_SUCCESS);
// Get adjoint of the slater-determinant
double det_inv_matrix_alpha[det_num_alpha][chbrclf_walk_num][chbrclf_elec_up_num][chbrclf_elec_up_num];
double det_inv_matrix_beta[det_num_beta][chbrclf_walk_num][chbrclf_elec_dn_num][chbrclf_elec_dn_num];
rc = qmckl_get_det_inv_matrix_alpha(context, &(det_inv_matrix_alpha[0][0][0][0]));
assert (rc == QMCKL_SUCCESS);
rc = qmckl_get_det_inv_matrix_beta(context, &(det_inv_matrix_beta[0][0][0][0]));
assert (rc == QMCKL_SUCCESS);
// Calculate the Kinetic energy
double kinetic_energy[chbrclf_walk_num];
rc = qmckl_get_kinetic_energy(context, &(kinetic_energy[0]));
assert (rc == QMCKL_SUCCESS);
#+end_src
** Potential energy
:PROPERTIES:
:Name: qmckl_compute_potential_energy
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
The potential energy is the sum of all the following terms
\[
PE = \mathcal{V}_{ee} + \mathcal{V}_{en} + \mathcal{V}_{nn}
\]
The potential for is calculated as the sum of single electron
contributions.
\[
\mathcal{V}_{ee} = \sum_{i=1}^{N_e}\sum_{j<i} \frac{1}{r_{ij}}
\]
\[
\mathcal{V}_{en} = - \sum_{i=1}^{N_e}\sum_{A=1}^{N_n}\frac{Z_A}{r_{iA}}
\]
\[
\mathcal{V}_{nn} = \sum_{A=1}^{N_n}\sum_{B<A}\frac{Z_A Z_B}{r_{AB}}
\]
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_potential_energy(qmckl_context context, double* const potential_energy);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_potential_energy(qmckl_context context, double * const potential_energy) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_exit_code rc;
if(!qmckl_electron_provided(context)) return QMCKL_NOT_PROVIDED;
if(!qmckl_nucleus_provided(context)) return QMCKL_NOT_PROVIDED;
rc = qmckl_provide_ao_basis_ao_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_mo_basis_mo_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_potential_energy(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->electron.walker.num * sizeof(double);
memcpy(potential_energy, ctx->local_energy.e_pot, sze);
return QMCKL_SUCCESS;
}
#+end_src
*** Provide
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_potential_energy(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_provide_potential_energy(qmckl_context context) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
qmckl_exit_code rc;
if(!(ctx->nucleus.provided)) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_electron",
NULL);
}
if(!(ctx->electron.provided)) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_electron",
NULL);
}
if (!ctx->ao_basis.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_ao_basis",
NULL);
}
if (!ctx->mo_basis.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_mo_basis",
NULL);
}
if (!ctx->det.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_mo_basis",
NULL);
}
rc = qmckl_provide_nucleus_repulsion(context);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_nucleus_repulsion",
NULL);
}
rc = qmckl_provide_ee_potential(context);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_ee_potential",
NULL);
}
rc = qmckl_provide_en_potential(context);
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_en_potential",
NULL);
}
/* Compute if necessary */
if (ctx->electron.walker.point.date > ctx->local_energy.e_pot_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->local_energy.e_pot);
ctx->local_energy.e_pot = NULL;
}
/* Allocate array */
if (ctx->local_energy.e_pot == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->electron.walker.num * sizeof(double);
double* e_pot = (double*) qmckl_malloc(context, mem_info);
if (e_pot == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_provide_potential_energy",
NULL);
}
ctx->local_energy.e_pot = e_pot;
}
if (ctx->det.type == 'G') {
rc = qmckl_compute_potential_energy(context,
ctx->electron.walker.num,
ctx->electron.num,
ctx->nucleus.num,
ctx->electron.ee_potential,
ctx->electron.en_potential,
ctx->nucleus.repulsion,
ctx->local_energy.e_pot);
} else {
return qmckl_failwith( context,
QMCKL_FAILURE,
"compute_potential_energy",
"Not yet implemented");
}
if (rc != QMCKL_SUCCESS) {
return rc;
}
ctx->local_energy.e_pot_date = ctx->date;
}
return QMCKL_SUCCESS;
}
#+end_src
*** Compute potential enregy
:PROPERTIES:
:Name: qmckl_compute_potential_energy
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
#+NAME: qmckl_compute_potential_energy_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 MOs |
| ~double~ | ~ee_potential[walk_num]~ | in | ee potential |
| ~double~ | ~en_potential[walk_num]~ | in | en potential |
| ~double~ | ~repulsion~ | in | en potential |
| ~double~ | ~e_pot[walk_num]~ | out | Potential energy |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_potential_energy_f(context, walk_num, &
elec_num, nucl_num, ee_potential, en_potential, repulsion, e_pot) &
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
double precision, intent(in) :: ee_potential(walk_num)
double precision, intent(in) :: en_potential(walk_num)
double precision, intent(in) :: repulsion
double precision, intent(inout) :: e_pot(walk_num)
integer*8 :: idet, iwalk, ielec, mo_id, imo
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
do iwalk = 1, walk_num
e_pot(iwalk) = ee_potential(iwalk) + en_potential(iwalk) + repulsion
end do
end function qmckl_compute_potential_energy_f
#+end_src
#+CALL: generate_c_header(table=qmckl_compute_potential_energy_args,rettyp=get_value("CRetType"),fname="qmckl_compute_potential_energy"))
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
qmckl_exit_code qmckl_compute_potential_energy (
const qmckl_context context,
const int64_t walk_num,
const int64_t elec_num,
const int64_t nucl_num,
const double* ee_potential,
const double* en_potential,
const double repulsion,
double* const e_pot );
#+end_src
#+CALL: generate_c_interface(table=qmckl_compute_potential_energy_args,rettyp=get_value("CRetType"),fname="qmckl_compute_potential_energy"))
#+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none
integer(c_int32_t) function qmckl_compute_potential_energy &
(context, walk_num, elec_num, nucl_num, ee_potential, en_potential, repulsion, e_pot) &
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
real (c_double ) , intent(in) :: ee_potential(walk_num)
real (c_double ) , intent(in) :: en_potential(walk_num)
real (c_double ) , intent(in) , value :: repulsion
real (c_double ) , intent(out) :: e_pot(walk_num)
integer(c_int32_t), external :: qmckl_compute_potential_energy_f
info = qmckl_compute_potential_energy_f &
(context, walk_num, elec_num, nucl_num, ee_potential, en_potential, repulsion, e_pot)
end function qmckl_compute_potential_energy
#+end_src
*** Test
#+begin_src c :tangle (eval c_test) :exports none
// Calculate the Potential energy
double potential_energy[chbrclf_walk_num];
rc = qmckl_get_potential_energy(context, &(potential_energy[0]));
assert (rc == QMCKL_SUCCESS);
#+end_src
** Local energy
:PROPERTIES:
:Name: qmckl_compute_local_energy
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
The local energy is the sum of kinetic and potential energies.
\[
E_L = KE + PE
\]
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_local_energy(qmckl_context context, double* const local_energy, const int64_t size_max);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_local_energy(qmckl_context context, double * const local_energy, const int64_t size_max) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_exit_code rc;
if(!qmckl_electron_provided(context)) return QMCKL_NOT_PROVIDED;
if(!qmckl_nucleus_provided(context)) return QMCKL_NOT_PROVIDED;
rc = qmckl_provide_ao_basis_ao_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_mo_basis_mo_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_local_energy(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
const int64_t sze = ctx->electron.walker.num;
if (size_max < sze) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_get_local_energy",
"input array too small");
}
memcpy(local_energy, ctx->local_energy.e_local, sze * sizeof(double));
return QMCKL_SUCCESS;
}
#+end_src
*** Provide
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_local_energy(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_provide_local_energy(qmckl_context context) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if(!(ctx->nucleus.provided)) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_electron",
NULL);
}
if(!(ctx->electron.provided)) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_electron",
NULL);
}
if (!ctx->ao_basis.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_ao_basis",
NULL);
}
if (!ctx->mo_basis.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_mo_basis",
NULL);
}
if (!ctx->det.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_mo_basis",
NULL);
}
qmckl_exit_code rc;
rc = qmckl_provide_kinetic_energy(context);
if(rc != QMCKL_SUCCESS){
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_kinetic_energy",
NULL);
}
rc = qmckl_provide_potential_energy(context);
if(rc != QMCKL_SUCCESS){
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_potential_energy",
NULL);
}
/* Compute if necessary */
if (ctx->electron.walker.point.date > ctx->local_energy.e_local_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->local_energy.e_local);
ctx->local_energy.e_local = NULL;
}
/* Allocate array */
if (ctx->local_energy.e_local == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->electron.walker.num * sizeof(double);
double* local_energy = (double*) qmckl_malloc(context, mem_info);
if (local_energy == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_local_energy",
NULL);
}
ctx->local_energy.e_local = local_energy;
}
if (ctx->det.type == 'G') {
rc = qmckl_compute_local_energy(context,
ctx->electron.walker.num,
ctx->local_energy.e_kin,
ctx->local_energy.e_pot,
ctx->local_energy.e_local);
} else {
return qmckl_failwith( context,
QMCKL_FAILURE,
"compute_local_energy",
"Not yet implemented");
}
if (rc != QMCKL_SUCCESS) {
return rc;
}
ctx->local_energy.e_local_date = ctx->date;
}
return QMCKL_SUCCESS;
}
#+end_src
*** Compute local enregy
:PROPERTIES:
:Name: qmckl_compute_local_energy
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
#+NAME: qmckl_compute_local_energy_args
| ~qmckl_context~ | ~context~ | in | Global state |
| ~int64_t~ | ~walk_num~ | in | Number of walkers |
| ~double~ | ~e_kin[walk_num]~ | in | e kinetic |
| ~double~ | ~e_pot[walk_num]~ | in | e potential |
| ~double~ | ~e_local[walk_num]~ | out | local energy |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_local_energy_f(context, walk_num, &
e_kin, e_pot, e_local) &
result(info)
use qmckl
implicit none
integer(qmckl_context) , intent(in) :: context
integer*8, intent(in) :: walk_num
double precision, intent(in) :: e_kin(walk_num)
double precision, intent(in) :: e_pot(walk_num)
double precision, intent(inout) :: e_local(walk_num)
integer*8 :: idet, iwalk, ielec, mo_id, imo
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
e_local = 0.0d0
do iwalk = 1, walk_num
e_local(iwalk) = e_local(iwalk) + e_kin(iwalk) + e_pot(iwalk)
end do
end function qmckl_compute_local_energy_f
#+end_src
#+CALL: generate_c_header(table=qmckl_compute_local_energy_args,rettyp=get_value("CRetType"),fname="qmckl_compute_local_energy"))
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
qmckl_exit_code qmckl_compute_local_energy (
const qmckl_context context,
const int64_t walk_num,
const double* e_kin,
const double* e_pot,
double* const e_local );
#+end_src
#+CALL: generate_c_interface(table=qmckl_compute_local_energy_args,rettyp=get_value("CRetType"),fname="qmckl_compute_local_energy"))
#+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none
integer(c_int32_t) function qmckl_compute_local_energy &
(context, walk_num, e_kin, e_pot, e_local) &
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
real (c_double ) , intent(in) :: e_kin(walk_num)
real (c_double ) , intent(in) :: e_pot(walk_num)
real (c_double ) , intent(out) :: e_local(walk_num)
integer(c_int32_t), external :: qmckl_compute_local_energy_f
info = qmckl_compute_local_energy_f &
(context, walk_num, e_kin, e_pot, e_local)
end function qmckl_compute_local_energy
#+end_src
*** Test
#+begin_src c :tangle (eval c_test) :exports none
// Calculate the Local energy
double local_energy[chbrclf_walk_num];
rc = qmckl_get_local_energy(context, &(local_energy[0]), chbrclf_walk_num);
assert (rc == QMCKL_SUCCESS);
#+end_src
** Drift vector
:PROPERTIES:
:Name: qmckl_compute_drift_vector
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
The drift vector is calculated as the ration of the gradient
with the determinant of the wavefunction.
\[
\mathbf{F} = 2 \frac{\nabla \Psi}{\Psi}
\]
*** Get
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_get_drift_vector(qmckl_context context, double* const drift_vector);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_get_drift_vector(qmckl_context context, double * const drift_vector) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_exit_code rc;
if(!qmckl_electron_provided(context)) return QMCKL_NOT_PROVIDED;
if(!qmckl_nucleus_provided(context)) return QMCKL_NOT_PROVIDED;
rc = qmckl_provide_ao_basis_ao_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_mo_basis_mo_vgl(context);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_provide_drift_vector(context);
if (rc != QMCKL_SUCCESS) return rc;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
size_t sze = ctx->electron.walker.num * ctx->electron.num * 3 * sizeof(double);
memcpy(drift_vector, ctx->local_energy.r_drift, sze);
return QMCKL_SUCCESS;
}
#+end_src
*** Provide
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
qmckl_exit_code qmckl_provide_drift_vector(qmckl_context context);
#+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_provide_drift_vector(qmckl_context context) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if(!(ctx->nucleus.provided)) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_electron",
NULL);
}
if(!(ctx->electron.provided)) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_electron",
NULL);
}
if (!ctx->ao_basis.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_ao_basis",
NULL);
}
if (!ctx->mo_basis.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_mo_basis",
NULL);
}
if (!ctx->det.provided) {
return qmckl_failwith( context,
QMCKL_NOT_PROVIDED,
"qmckl_mo_basis",
NULL);
}
/* Compute if necessary */
if (ctx->electron.walker.point.date > ctx->local_energy.r_drift_date) {
if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
free(ctx->local_energy.r_drift);
ctx->local_energy.r_drift = NULL;
}
/* Allocate array */
if (ctx->local_energy.r_drift == NULL) {
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = ctx->electron.walker.num * ctx->electron.num * 3 * sizeof(double);
double* r_drift = (double*) qmckl_malloc(context, mem_info);
if (r_drift == NULL) {
return qmckl_failwith( context,
QMCKL_ALLOCATION_FAILED,
"qmckl_r_drift",
NULL);
}
ctx->local_energy.r_drift = r_drift;
}
qmckl_exit_code rc;
if (ctx->det.type == 'G') {
rc = qmckl_compute_drift_vector(context,
ctx->electron.walker.num,
ctx->det.det_num_alpha,
ctx->det.det_num_beta,
ctx->electron.up_num,
ctx->electron.down_num,
ctx->electron.num,
ctx->det.mo_index_alpha,
ctx->det.mo_index_beta,
ctx->mo_basis.mo_num,
ctx->mo_basis.mo_vgl,
ctx->det.det_inv_matrix_alpha,
ctx->det.det_inv_matrix_beta,
ctx->local_energy.r_drift);
} else {
return qmckl_failwith( context,
QMCKL_FAILURE,
"compute_drift_vector",
"Not yet implemented");
}
if (rc != QMCKL_SUCCESS) {
return rc;
}
ctx->local_energy.r_drift_date = ctx->date;
}
return QMCKL_SUCCESS;
}
#+end_src
*** Compute drift vector
:PROPERTIES:
:Name: qmckl_compute_drift_vector
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
#+NAME: qmckl_compute_drift_vector_args
| ~qmckl_context~ | ~context~ | in | Global state |
| ~int64_t~ | ~walk_num~ | in | Number of walkers |
| ~int64_t~ | ~det_num_alpha~ | in | Number of determinants |
| ~int64_t~ | ~det_num_beta~ | in | Number of determinants |
| ~int64_t~ | ~alpha_num~ | in | Number of electrons |
| ~int64_t~ | ~beta_num~ | in | Number of electrons |
| ~int64_t~ | ~elec_num~ | in | Number of electrons |
| ~int64_t~ | ~mo_index_alpha[det_num_alpha][walk_num][alpha_num]~ | in | MO indices for electrons |
| ~int64_t~ | ~mo_index_beta[det_num_beta][walk_num][beta_num]~ | in | MO indices for electrons |
| ~int64_t~ | ~mo_num~ | in | Number of MOs |
| ~double~ | ~mo_vgl[5][elec_num][mo_num]~ | in | Value, gradients and Laplacian of the MOs |
| ~double~ | ~det_inv_matrix_alpha[det_num_alpha][walk_num][alpha_num][alpha_num]~ | in | Value, gradients and Laplacian of the Det |
| ~double~ | ~det_inv_matrix_beta[det_num_beta][walk_num][beta_num][beta_num]~ | in | Value, gradients and Laplacian of the Det |
| ~double~ | ~r_drift[walk_num][elec_num][3]~ | out | Kinetic energy |
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_drift_vector_f(context, walk_num, &
det_num_alpha, det_num_beta, alpha_num, beta_num, elec_num, mo_index_alpha, mo_index_beta, &
mo_num, mo_vgl, det_inv_matrix_alpha, det_inv_matrix_beta, r_drift) &
result(info)
use qmckl
implicit none
integer(qmckl_context) , intent(in) :: context
integer*8, intent(in) :: walk_num
integer*8, intent(in) :: det_num_alpha
integer*8, intent(in) :: det_num_beta
integer*8, intent(in) :: alpha_num
integer*8, intent(in) :: beta_num
integer*8, intent(in) :: elec_num
integer*8, intent(in) :: mo_num
integer*8, intent(in) :: mo_index_alpha(alpha_num, walk_num, det_num_alpha)
integer*8, intent(in) :: mo_index_beta(beta_num, walk_num, det_num_beta)
double precision, intent(in) :: mo_vgl(mo_num, elec_num, 5)
double precision, intent(in) :: det_inv_matrix_alpha(alpha_num, alpha_num, walk_num, det_num_alpha)
double precision, intent(in) :: det_inv_matrix_beta(beta_num, beta_num, walk_num, det_num_beta)
double precision, intent(inout) :: r_drift(3,elec_num,walk_num)
integer*8 :: idet, iwalk, ielec, mo_id, imo
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 (alpha_num <= 0) then
info = QMCKL_INVALID_ARG_3
return
endif
if (beta_num < 0) then
info = QMCKL_INVALID_ARG_4
return
endif
if (elec_num <= 0) then
info = QMCKL_INVALID_ARG_5
return
endif
r_drift = 0.0d0
do idet = 1, det_num_alpha
do iwalk = 1, walk_num
! Alpha part
do imo = 1, alpha_num
do ielec = 1, alpha_num
mo_id = mo_index_alpha(imo, iwalk, idet)
r_drift(1,ielec,iwalk) = r_drift(1,ielec,iwalk) + 2.0d0 * det_inv_matrix_alpha(imo, ielec, iwalk, idet) * &
mo_vgl(mo_id, ielec, 2)
r_drift(2,ielec,iwalk) = r_drift(2,ielec,iwalk) + 2.0d0 * det_inv_matrix_alpha(imo, ielec, iwalk, idet) * &
mo_vgl(mo_id, ielec, 3)
r_drift(3,ielec,iwalk) = r_drift(3,ielec,iwalk) + 2.0d0 * det_inv_matrix_alpha(imo, ielec, iwalk, idet) * &
mo_vgl(mo_id, ielec, 4)
end do
end do
! Beta part
do imo = 1, beta_num
do ielec = 1, beta_num
mo_id = mo_index_beta(imo, iwalk, idet)
r_drift(1,alpha_num + ielec,iwalk) = r_drift(1,alpha_num + ielec,iwalk) + &
2.0d0 * det_inv_matrix_beta(imo, ielec, iwalk, idet) * &
mo_vgl(mo_id, alpha_num + ielec, 2)
r_drift(2,alpha_num + ielec,iwalk) = r_drift(2,alpha_num + ielec,iwalk) + &
2.0d0 * det_inv_matrix_beta(imo, ielec, iwalk, idet) * &
mo_vgl(mo_id, alpha_num + ielec, 3)
r_drift(3,alpha_num + ielec,iwalk) = r_drift(3,alpha_num + ielec,iwalk) + &
2.0d0 * det_inv_matrix_beta(imo, ielec, iwalk, idet) * &
mo_vgl(mo_id, alpha_num + ielec, 4)
end do
end do
end do
end do
end function qmckl_compute_drift_vector_f
#+end_src
#+CALL: generate_c_header(table=qmckl_compute_drift_vector_args,rettyp=get_value("CRetType"),fname="qmckl_compute_drift_vector"))
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
qmckl_exit_code qmckl_compute_drift_vector (
const qmckl_context context,
const int64_t walk_num,
const int64_t det_num_alpha,
const int64_t det_num_beta,
const int64_t alpha_num,
const int64_t beta_num,
const int64_t elec_num,
const int64_t* mo_index_alpha,
const int64_t* mo_index_beta,
const int64_t mo_num,
const double* mo_vgl,
const double* det_inv_matrix_alpha,
const double* det_inv_matrix_beta,
double* const r_drift );
#+end_src
#+CALL: generate_c_interface(table=qmckl_compute_drift_vector_args,rettyp=get_value("CRetType"),fname="qmckl_compute_drift_vector"))
#+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none
integer(c_int32_t) function qmckl_compute_drift_vector &
(context, &
walk_num, &
det_num_alpha, &
det_num_beta, &
alpha_num, &
beta_num, &
elec_num, &
mo_index_alpha, &
mo_index_beta, &
mo_num, &
mo_vgl, &
det_inv_matrix_alpha, &
det_inv_matrix_beta, &
r_drift) &
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 :: det_num_alpha
integer (c_int64_t) , intent(in) , value :: det_num_beta
integer (c_int64_t) , intent(in) , value :: alpha_num
integer (c_int64_t) , intent(in) , value :: beta_num
integer (c_int64_t) , intent(in) , value :: elec_num
integer (c_int64_t) , intent(in) :: mo_index_alpha(alpha_num,walk_num,det_num_alpha)
integer (c_int64_t) , intent(in) :: mo_index_beta(beta_num,walk_num,det_num_beta)
integer (c_int64_t) , intent(in) , value :: mo_num
real (c_double ) , intent(in) :: mo_vgl(mo_num,elec_num,5)
real (c_double ) , intent(in) :: det_inv_matrix_alpha(alpha_num,alpha_num,walk_num,det_num_alpha)
real (c_double ) , intent(in) :: det_inv_matrix_beta(beta_num,beta_num,walk_num,det_num_beta)
real (c_double ) , intent(out) :: r_drift(3,elec_num,walk_num)
integer(c_int32_t), external :: qmckl_compute_drift_vector_f
info = qmckl_compute_drift_vector_f &
(context, &
walk_num, &
det_num_alpha, &
det_num_beta, &
alpha_num, &
beta_num, &
elec_num, &
mo_index_alpha, &
mo_index_beta, &
mo_num, &
mo_vgl, &
det_inv_matrix_alpha, &
det_inv_matrix_beta, &
r_drift)
end function qmckl_compute_drift_vector
#+end_src
*** Test
#+begin_src c :tangle (eval c_test) :exports none
// Calculate the Drift vector
double drift_vector[chbrclf_walk_num][chbrclf_elec_num][3];
rc = qmckl_get_drift_vector(context, &(drift_vector[0][0][0]));
assert (rc == QMCKL_SUCCESS);
#+end_src
* End of files :noexport:
#+begin_src c :tangle (eval h_private_type)
#endif
#+end_src
*** Test
#+begin_src c :tangle (eval c_test)
rc = qmckl_context_destroy(context);
assert (rc == QMCKL_SUCCESS);
return 0;
}
#+end_src
*** Compute file names
#+begin_src emacs-lisp
; The following is required to compute the file names
(setq pwd (file-name-directory buffer-file-name))
(setq name (file-name-nondirectory (substring buffer-file-name 0 -4)))
(setq f (concat pwd name "_f.f90"))
(setq fh (concat pwd name "_fh.f90"))
(setq c (concat pwd name ".c"))
(setq h (concat name ".h"))
(setq h_private (concat name "_private.h"))
(setq c_test (concat pwd "test_" name ".c"))
(setq f_test (concat pwd "test_" name "_f.f90"))
; Minted
(require 'ox-latex)
(setq org-latex-listings 'minted)
(add-to-list 'org-latex-packages-alist '("" "listings"))
(add-to-list 'org-latex-packages-alist '("" "color"))
#+end_src
# -*- mode: org -*-
# vim: syntax=c
Reference in New Issue View Git Blame Copy Permalink