mirror of
https://github.com/TREX-CoE/qmckl.git
synced 2024-11-05 05:34:06 +01:00
58 KiB
58 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:
- Drift Vector - \[F(x)\]
\[ F(x) = 2\frac{\nabla \Psi(r)}{\Psi(r)} \]
- 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.
- 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][elec_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;
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.
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);Computation
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) \]
Get
qmckl_exit_code qmckl_get_kinetic_energy(qmckl_context context, double* const kinetic_energy);Provide
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_fqmckl_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 );Test
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<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
qmckl_exit_code qmckl_get_potential_energy(qmckl_context context, double* const potential_energy);Provide
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_fqmckl_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 );Test
Local energy
The local energy is the sum of kinetic and potential energies.
\[ E_L = KE + PE \]
Get
qmckl_exit_code qmckl_get_local_energy(qmckl_context context, double* const local_energy, const int64_t size_max);Provide
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_fqmckl_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 );Test
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} \]
Get
qmckl_exit_code qmckl_get_drift_vector(qmckl_context context, double* const drift_vector);Provide
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_fqmckl_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 );