QMC/Diffusion.cpp at master · jorgehog/QMC · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
/*
 * File:   Diffusion.cpp
 * Author: jorgehog
 *
 * Created on 16. april 2012, 14:03
 */

#include "QMCheaders.h"

Diffusion::Diffusion(int n_p, int dim, double timestep, long random_seed, double D) {
    this->n_p = n_p;
    this->dim = dim;
    this->timestep = timestep;
    this->random_seed = random_seed;
    this->D = D;
    this->std = sqrt(2 * D * timestep);
}

double Diffusion::call_RNG() {

    return ran3(&random_seed);
}

double Diffusion::get_new_pos(const Walker* walker, int i, int j) {
    return gaussian_deviate(&random_seed) * std;
}

Simple::Simple(int n_p, int dim, double timestep, long random_seed, double D)
: Diffusion(n_p, dim, timestep, random_seed, D) {

}

double Simple::get_new_pos(const Walker* walker, int i, int j) {
    return Diffusion::get_new_pos(walker, i, j);
}

double Simple::get_g_ratio(const Walker* walker_post, const Walker* walker_pre, int particle) const {
    return 1;
}

double Simple::get_GBfunc(Walker* walker_pre, Walker* walker_post, double E_T) const {
    double Vx = qmc->get_system_ptr()->get_potential_energy(walker_pre);
    double Vy = qmc->get_system_ptr()->get_potential_energy(walker_post);

    return exp(-(0.5 * (Vy + Vx) - E_T) * timestep);


}

Fokker_Planck::Fokker_Planck(int n_p, int dim, double timestep, long random_seed, double D)
: Diffusion(n_p, dim, timestep, random_seed, D) {

}

double Fokker_Planck::get_new_pos(const Walker* walker, int i, int j) {
    return D * timestep * walker->qforce(i, j) + Diffusion::get_new_pos(walker, i, j);
}

double Fokker_Planck::get_g_ratio(const Walker* walker_post, const Walker* walker_pre, int particle) const {

    double g_ratio = 0;
    for (int j = 0; j < dim; j++) {
        g_ratio += 0.5 * (walker_pre->qforce(particle, j) + walker_post->qforce(particle, j))*
                (D * timestep * 0.5 * (walker_pre->qforce(particle, j) - walker_post->qforce(particle, j))
                - walker_post->r(particle, j) + walker_pre->r(particle, j));
    }

    return exp(g_ratio);
}

double Fokker_Planck::get_GBfunc(Walker* walker_pre, Walker* walker_post, double E_T) const {

    double Ex = qmc->calculate_local_energy(walker_pre);
    double Ey = qmc->calculate_local_energy(walker_post);

    return exp(-(0.5*(Ey + Ex) - E_T)*timestep*qmc->get_accepted_ratio());
}