484 14 Quantum Monte Carlo Methods
// update energies
energy += delta_e;
energy2 += delta_e*delta_e;
}
} // end of loop over MC trials
cout <<"variational parameter= "<< alpha
<<" accepted steps= "<< accept << endl;
// update the energy average and its squared
cumulative_e[variate] = energy/number_cycles;
cumulative_e2[variate] = energy2/number_cycles;
} // end of loop over variational steps
free_matrix((void) r_old);// free memory
free_matrix((void) r_new);// free memory
} // end mc_sampling function
The wave function is in turn defined in the next function. Herewe limit ourselves to a function
which consists only of the product of single-particle wave functions.
// Function to compute the squared wave function, simplest form
doublewave_function(double *r,doublealpha,intdimension,intnumber_particles)
{
inti, j, k;
doublewf, argument, r_single_particle, r_12;
argument = wf = 0;
for(i = 0; i < number_particles; i++){
r_single_particle = 0;
for(j = 0; j < dimension; j++){
r_single_particle += r[i][j]r[i][j];
}
argument += sqrt(r_single_particle);
}
wf = exp(-argument*alpha) ;
returnwf;
}
Finally, the local energy is computed using a numerical derivation for the kinetic energy. We
use the familiar expression derived in Eq. (3.4), that is
f 0 ′′=
fh− 2 f 0 +f−h
h^2 ,in order to compute
−1
2 ψT(R)
∇^2 ψT(R). (14.23)The variablehis a chosen step length. For helium, since it is rather easy toevaluate the
local energy, the above is an unnecessary complication. However, for many-electron or other
many-particle systems, the derivation of a closed-form expression for the kinetic energy can
be quite involved, and the numerical evaluation of the kinetic energy using Eq. (3.4) may
result in a simpler code and/or even a faster one.
// Function to calculate the local energy with num derivative
doublelocal_energy(double**r,doublealpha,doublewfold,intdimension,
intnumber_particles,intcharge)
{
inti, j , k;
doublee_local, wfminus, wfplus, e_kinetic, e_potential, r_12,