UP | HOME

Local Energy

Table of Contents

1 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

1.1 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;
  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;

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.

2 Computation

2.1 Kinetic energy

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

2.1.1 Get

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

2.1.2 Provide

2.1.3 Compute kinetic enregy

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

qmckl_exit_code qmckl_compute_kinetic_energy (
      const context qmckl_context,
      const walk_num int64_t,
      const det_num_alpha int64_t,
      const det_num_beta int64_t,
      const alpha_num int64_t,
      const beta_num int64_t,
      const elec_num int64_t,
      const mo_index_alpha* int64_t,
      const mo_index_beta* int64_t,
      const mo_num int64_t,
      const mo_vgl* double,
      const det_value_alpha* double,
      const det_value_beta* double,
      const det_inv_matrix_alpha* double,
      const det_inv_matrix_beta* double,
      e_kin* const double );

2.1.4 Test

2.2 Potential energy

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

\[ \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

2.2.1 Get

qmckl_exit_code qmckl_get_potential_energy(qmckl_context context, double* const potential_energy);

2.2.2 Provide

2.2.3 Compute potential enregy

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

qmckl_exit_code qmckl_compute_potential_energy (
      const context qmckl_context,
      const walk_num int64_t,
      const elec_num int64_t,
      const nucl_num int64_t,
      const ee_potential* double,
      const en_potential* double,
      const repulsion double,
      e_pot* const double );

2.2.4 Test

2.3 Local energy

The local energy is the sum of kinetic and potential energies.

\[ E_L = KE + PE \]

2.3.1 Get

qmckl_exit_code qmckl_get_local_energy(qmckl_context context, double* const local_energy, const int64_t size_max);

2.3.2 Provide

2.3.3 Compute local enregy

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

qmckl_exit_code qmckl_compute_local_energy (
      const context qmckl_context,
      const walk_num int64_t,
      const e_kin* double,
      const e_pot* double,
      e_local* const double );

2.3.4 Test

2.4 Drift vector

The drift vector is calculated as the ration of the gradient with the determinant of the wavefunction.

\[ \mathbf{F} = 2 \frac{\nabla \Psi}{\Psi} \]

2.4.1 Get

qmckl_exit_code qmckl_get_drift_vector(qmckl_context context, double* const drift_vector);

2.4.2 Provide

2.4.3 Compute drift vector

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

qmckl_exit_code qmckl_compute_drift_vector (
      const context qmckl_context,
      const walk_num int64_t,
      const det_num_alpha int64_t,
      const det_num_beta int64_t,
      const alpha_num int64_t,
      const beta_num int64_t,
      const elec_num int64_t,
      const mo_index_alpha* int64_t,
      const mo_index_beta* int64_t,
      const mo_num int64_t,
      const mo_vgl* double,
      const det_inv_matrix_alpha* double,
      const det_inv_matrix_beta* double,
      r_drift* const double );

2.4.4 Test

Author: TREX CoE

Created: 2022-08-08 Mon 16:19

Validate