mirror of
https://github.com/TREX-CoE/qmckl.git
synced 2024-12-23 12:57:20 +01:00
1739 lines
58 KiB
Org Mode
1739 lines
58 KiB
Org Mode
#+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);
|
|
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
|