Skip to content

Commit 7f2f126

Browse files
yaozhenghangmayaozhenghangma
committed
Ising model
1 parent 3e39a6c commit 7f2f126

File tree

5 files changed

+253
-4
lines changed

5 files changed

+253
-4
lines changed

‎include/Ising.hpp‎

Lines changed: 227 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,227 @@
1+
#ifndef HYBRID_PROGRAMMING_ISING_HPP
2+
#define HYBRID_PROGRAMMING_ISING_HPP
3+
#include <iostream>
4+
#include <numeric>
5+
#include <random>
6+
#include <stdexcept>
7+
#include <vector>
8+
9+
template <class T>
10+
using Vector2d = std::vector<std::vector<T>>;
11+
12+
const double one_over_KB = 11.60450501;
13+
14+
class Lattice2d{
15+
public:
16+
Vector2d<int> spin;
17+
Lattice2d(){}
18+
Lattice2d(int nx, int ny) {
19+
auto vector1d = std::vector<int>(nx, 1);
20+
this->spin = Vector2d<int>(ny, vector1d);
21+
this->max_index_x_ = nx-1;
22+
this->max_index_y_ = ny-1;
23+
}
24+
25+
std::tuple<int, int, int, int> neighbors(int ix, int iy) {
26+
if(ix == 0) {
27+
if(iy == 0) {
28+
return {this->max_index_x_, ix+1, this->max_index_y_, iy+1};
29+
} else if(iy == this->max_index_y_) {
30+
return {this->max_index_x_, ix+1, iy-1, 0};
31+
} else {
32+
return {this->max_index_x_, ix+1, iy-1, iy+1};
33+
}
34+
} else if(ix == this->max_index_x_) {
35+
if(iy == 0) {
36+
return {ix-1, 0, this->max_index_y_, iy+1};
37+
} else if(iy == this->max_index_y_) {
38+
return {ix-1, 0, iy-1, 0};
39+
} else {
40+
return {ix-1, 0, iy-1, iy+1};
41+
}
42+
} else {
43+
if(iy == 0) {
44+
return {ix-1, ix+1, this->max_index_y_, iy+1};
45+
} else if(iy == this->max_index_y_) {
46+
return {ix-1, ix+1, iy-1, 0};
47+
} else {
48+
return {ix-1, ix+1, iy-1, iy+1};
49+
}
50+
}
51+
}
52+
53+
private:
54+
int max_index_x_ = -1;
55+
int max_index_y_ = -1;
56+
};
57+
58+
class MonteCarlo {
59+
public:
60+
int num_flip = 0;
61+
int num_relax = 0;
62+
int num_count = 0;
63+
Lattice2d lattice;
64+
double energy = 0;
65+
double Cv = 0;
66+
67+
MonteCarlo(double J, int nx, int ny) {
68+
this->nx_ = nx;
69+
this->ny_ = ny;
70+
this->lattice = Lattice2d(nx, ny);
71+
this->SetJ(J);
72+
}
73+
74+
double GetJ();
75+
double GetT();
76+
int SetJ(double J);
77+
int SetT(double T);
78+
int Run();
79+
80+
private:
81+
double J_ = 0;
82+
double T_ = 0; //temperature
83+
double current_energy_ = 0;
84+
double total_energy_ = 0;
85+
double total_energy_square_ = 0;
86+
int nx_ = 0;
87+
int ny_ = 0;
88+
double kEnergy_positive_4_ = 0;
89+
double kEnergy_positive_2_ = 0;
90+
double kEnergy_0_ = 0;
91+
double kEnergy_negative_2_ = 0;
92+
double kEnergy_negative_4_ = 0;
93+
double kP_positive_4_ = 0;
94+
double kP_positive_2_ = 0;
95+
double kP_0_ = 0;
96+
double kP_negative_2_ = 0;
97+
double kP_negative_4_ = 0;
98+
99+
double CalculateEnergy();
100+
int MonteCarloStep();
101+
int Flip(int x, int y);
102+
std::tuple<double, double> EnergyAndProbability(int spin_sum);
103+
std::tuple<int, int> RandomSite();
104+
double RandomFloat();
105+
};
106+
107+
double MonteCarlo::GetJ() {return this->J_;}
108+
double MonteCarlo::GetT() {return this->T_;}
109+
110+
int MonteCarlo::SetJ(double J) {
111+
this->J_ = J;
112+
this->kEnergy_positive_4_ = -8*this->J_;
113+
this->kEnergy_positive_2_ = -4*this->J_;
114+
this->kEnergy_0_ = 0;
115+
this->kEnergy_negative_2_ = 4*this->J_;
116+
this->kEnergy_negative_4_ = 8*this->J_;
117+
this->current_energy_ = CalculateEnergy();
118+
this->SetT(this->T_);
119+
return 0;
120+
}
121+
122+
int MonteCarlo::SetT(double T) {
123+
this->T_ = T;
124+
double one_over_T = 1.0/T;
125+
this->kP_positive_4_ = exp(-this->kEnergy_positive_4_ * one_over_KB * one_over_T);
126+
this->kP_positive_2_ = exp(-this->kEnergy_positive_2_ * one_over_KB * one_over_T);
127+
this->kP_0_ = exp(-this->kEnergy_0_ * one_over_KB * one_over_T);
128+
this->kP_negative_2_ = exp(-this->kEnergy_negative_2_ * one_over_KB * one_over_T);
129+
this->kP_negative_4_ = exp(-this->kEnergy_negative_4_ * one_over_KB * one_over_T);
130+
return 0;
131+
}
132+
133+
int MonteCarlo::Run() {
134+
for(int i=0; i<this->num_relax; i++) {
135+
this->MonteCarloStep();
136+
}
137+
for(int i=0; i<this->num_count; i++) {
138+
this->MonteCarloStep();
139+
this->total_energy_ += this->current_energy_;
140+
this->total_energy_square_ += this->current_energy_ * this->current_energy_;
141+
}
142+
143+
this->energy = this->total_energy_ / (this->nx_ * this->ny_ * this->num_count);
144+
this->Cv = (this->total_energy_square_ * this->num_count - this->total_energy_ * this->total_energy_) *
145+
one_over_KB /
146+
(this->nx_ * this->ny_ * this->T_ * this->T_ * this->num_count * this->num_count);
147+
return 0;
148+
}
149+
150+
double MonteCarlo::CalculateEnergy() {
151+
this->current_energy_ = 0;
152+
int x_previous, x_next, y_previous, y_next;
153+
for(int ix=0; ix<this->nx_; ix++) {
154+
for(int iy=0; iy<this->ny_; iy++) {
155+
std::tie(x_previous, x_next, y_previous, y_next) = this->lattice.neighbors(ix, iy);
156+
this->current_energy_ += 0.5 * this->J_ * this->lattice.spin[ix][iy] * this->lattice.spin[x_previous][iy];
157+
this->current_energy_ += 0.5 * this->J_ * this->lattice.spin[ix][iy] * this->lattice.spin[x_next][iy];
158+
this->current_energy_ += 0.5 * this->J_ * this->lattice.spin[ix][iy] * this->lattice.spin[ix][y_previous];
159+
this->current_energy_ += 0.5 * this->J_ * this->lattice.spin[ix][iy] * this->lattice.spin[ix][y_next];
160+
}
161+
}
162+
return this->current_energy_;
163+
}
164+
165+
int MonteCarlo::MonteCarloStep() {
166+
int x, y;
167+
for(int i=0; i<this->num_flip; i++) {
168+
std::tie(x, y) = RandomSite();
169+
this->Flip(x, y);
170+
}
171+
return 0;
172+
}
173+
174+
std::tuple<double, double> MonteCarlo::EnergyAndProbability(int spin_sum) {
175+
if(spin_sum == 4) {
176+
return {this->kEnergy_positive_4_, this->kP_positive_4_};
177+
} else if(spin_sum == 2) {
178+
return {this->kEnergy_positive_2_, this->kP_positive_2_};
179+
} else if(spin_sum == 0) {
180+
return {this->kEnergy_0_, this->kP_0_};
181+
} else if(spin_sum == -2) {
182+
return {this->kEnergy_negative_2_, this->kP_negative_2_};
183+
} else if(spin_sum == -4) {
184+
return {this->kEnergy_negative_4_, this->kP_negative_4_};
185+
} else {
186+
throw std::out_of_range("Unknown value of spin summation.");
187+
}
188+
}
189+
190+
int MonteCarlo::Flip(int x, int y) {
191+
auto[x_previous, x_next, y_previous, y_next] = this->lattice.neighbors(x, y);
192+
int spin_sum =
193+
this->lattice.spin[x_previous][y] +
194+
this->lattice.spin[x_next][y] +
195+
this->lattice.spin[x][y_previous] +
196+
this->lattice.spin[x][y_next];
197+
spin_sum *= this->lattice.spin[x][y];
198+
199+
try {
200+
auto[delta_energy, probability] = this->EnergyAndProbability(spin_sum);
201+
if( this->RandomFloat() < probability) {
202+
this->lattice.spin[x][y] *= -1;
203+
this->current_energy_ += delta_energy;
204+
}
205+
} catch(std::out_of_range & err) {
206+
std::cerr << err.what() << std::endl;
207+
exit(-1);
208+
}
209+
return 0;
210+
}
211+
212+
std::tuple<int, int> MonteCarlo::RandomSite() {
213+
static std::random_device rd;
214+
static std::mt19937 engine(rd());
215+
static std::uniform_int_distribution<int> int_distribution_x(0, this->nx_-1);
216+
static std::uniform_int_distribution<int> int_distribution_y(0, this->ny_-1);
217+
return {int_distribution_x(engine), int_distribution_y(engine)};
218+
}
219+
220+
double MonteCarlo::RandomFloat() {
221+
static std::random_device rd;
222+
static std::mt19937 engine(rd());
223+
static std::uniform_real_distribution<double> float_distribution(0, 1);
224+
return float_distribution(engine);
225+
}
226+
227+
#endif //HYBRID_PROGRAMMING_ISING_HPP

‎setup.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -14,7 +14,7 @@
1414

1515
setup(
1616
name="hybrid_example",
17-
version="0.2.2",
17+
version="0.3.0",
1818
description="Tutorials about hybrid programming using C++ and Python",
1919
author="Yaozhenghang Ma",
2020
author_email="Yaozhenghang.Ma@gmail.com",

‎src/hybrid_example/CMakeLists.txt‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -12,5 +12,7 @@ find_package(pybind11 CONFIG REQUIRED PATHS ${pybind11_dir})
1212

1313
# target
1414
pybind11_add_module(example_cpp_lib example_export.cpp)
15+
pybind11_add_module(Ising_lib Ising_export.cpp)
1516

16-
install(TARGETS example_cpp_lib DESTINATION .)
17+
install(TARGETS example_cpp_lib DESTINATION .)
18+
install(TARGETS Ising_lib DESTINATION .)

‎src/hybrid_example/Ising_export.cpp‎

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,19 @@
1+
#include <pybind11/pybind11.h>
2+
#include <pybind11/stl.h>
3+
4+
#include "Ising.hpp"
5+
6+
namespace py = pybind11;
7+
8+
PYBIND11_MODULE(Ising_lib, m) {
9+
py::class_<MonteCarlo>(m, "MonteCarlo")
10+
.def_readwrite("num_flip", &MonteCarlo::num_flip)
11+
.def_readwrite("num_relax", &MonteCarlo::num_relax)
12+
.def_readwrite("num_count", &MonteCarlo::num_count)
13+
.def_readwrite("energy", &MonteCarlo::energy)
14+
.def_readwrite("Cv", &MonteCarlo::Cv)
15+
.def_property("J", &MonteCarlo::GetJ, &MonteCarlo::SetJ)
16+
.def_property("T", &MonteCarlo::GetT, &MonteCarlo::SetT)
17+
.def(py::init<double, int, int>())
18+
.def("Run", &MonteCarlo::Run);
19+
}

‎src/hybrid_example/__init__.py‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,5 @@
11
__author__ = 'Yaozhenghang Ma'
2-
__version__ = '0.0.1'
2+
__version__ = '0.3.0'
33

4-
from hybrid_example.example_cpp_lib import c_vector
4+
from hybrid_example.example_cpp_lib import c_vector
5+
from hybrid_example.Ising_lib import MonteCarlo

0 commit comments

Comments
 (0)