Local Energy

1. Context
- 1.1. Data structure
- 1.2. Access functions
2. Computation

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.

1.2 Access functions

When all the data for the local energy have been provided, the following function returns true.

bool      qmckl_local_energy_provided           (const qmckl_context context);

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