14.5 Variational Monte Carlo for atoms 483
The implementation of the brute force Metropolis algorithmis shown in the next function.
Here we have a loop over the variational variablesα. It calls two functions, one to compute
the wave function and one to update the local energy.
// Monte Carlo sampling with the Metropolis algorithm
voidmc_sampling(intdimension,intnumber_particles,intcharge,
intmax_variations,
intthermalization,intnumber_cycles,doublestep_length,
doublecumulative_e,doublecumulative_e2)
{
intcycles, variate, accept, dim, i, j;
longidum;
doublewfnew, wfold, alpha, energy, energy2, delta_e;
doubler_old,r_new;
alpha = 0.5charge;
idum=-1;
// allocate matrices which contain the position of the particles
r_old = (double) matrix( number_particles, dimension,sizeof(double));
r_new = (double) matrix( number_particles, dimension,sizeof(double));
for(i = 0; i < number_particles; i++){
for( j=0; j < dimension; j++){
r_old[i][j] = r_new[i][j] = 0;
}
}
// loop over variational parameters
for(variate=1; variate <= max_variations; variate++){
// initialisations of variational parameters and energies
alpha += 0.1;
energy = energy2 = 0; accept =0; delta_e=0;
// initial trial position, note calling with alpha
// and in three dimensions
for(i = 0; i < number_particles; i++){
for( j=0; j < dimension; j++){
r_old[i][j] = step_length(ran1(&idum)-0.5);
}
}
wfold = wave_function(r_old, alpha, dimension, number_particles);
// loop over monte carlo cycles
for(cycles = 1; cycles <= number_cycles+thermalization; cycles++){
// new position
for(i = 0; i < number_particles; i++){
for( j=0; j < dimension; j++){
r_new[i][j] = r_old[i][j]+step_length(ran1(&idum)-0.5);
}
}
wfnew = wave_function(r_new, alpha, dimension, number_particles);
// Metropolis test
if(ran1(&idum) <= wfnewwfnew/wfold/wfold ){
for(i = 0; i < number_particles; i++){
for( j=0; j < dimension; j++){
r_old[i][j]=r_new[i][j];
}
}
wfold = wfnew;
accept = accept+1;
}
// compute local energy
if( cycles > thermalization ){
delta_e = local_energy(r_old, alpha, wfold, dimension,
number_particles, charge);