1
0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-11-19 20:42:50 +01:00
qmckl/org/qmckl_local_energy.org
2021-10-11 16:12:14 +02:00

17 KiB

Local Energy

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)} \]

  1. 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.

  1. 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.

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] The drift vector
y_move [3][walk_num] The diffusion move
accep_prob [walk_num] The acceptance probability

Data structure

typedef struct qmckl_local_energy_struct {
double  * e_kin;
double  * e_pot;
double  * e_local;
double  * accep_prob;
double  * r_drift;
double  * y_move;
int64_t   e_kin_date;
int64_t   e_pot_date;
int64_t   e_local_date;
int64_t   accep_prob_date;
int64_t   r_drift_date;
int64_t   y_move_date;

int32_t   uninitialized;
bool      provided;
} qmckl_local_energy_struct;

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.

Computation

Kinetic energy

Get

qmckl_exit_code qmckl_get_kinetic_energy(qmckl_context context, double* const kinetic_energy);

Provide

Compute alpha

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][walk_num][elec_num][mo_num] in Value, gradients and Laplacian of the MOs
double det_vgl_alpha[det_num_alpha][walk_num][5][alpha_num][alpha_num] in Value, gradients and Laplacian of the Det
double det_vgl_beta[det_num_beta][walk_num][5][beta_num][beta_num] in Value, gradients and Laplacian of the Det
double e_kin[walk_num] out Kinetic energy
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_vgl_alpha, det_vgl_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, walk_num, 5)
double precision, intent(in)      :: det_vgl_alpha(alpha_num, alpha_num, 5, walk_num, det_num_alpha)
double precision, intent(in)      :: det_vgl_beta(beta_num, beta_num, 5, walk_num, det_num_beta)
double precision, intent(inout)   :: e_kin(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

idet = 1
e_kin = 0.0d0
do iwalk = 1, walk_num
! Alpha part
do ielec = 1, alpha_num
  mo_id = mo_index_beta(ielec, iwalk, idet)
end do
end do

end function qmckl_compute_kinetic_energy_f
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_vgl_alpha,
const double* det_vgl_beta,
double* const e_kin );

Test