问题
我正在尝试对 2D Ising 模型进行大都市模拟。
基本上，对于每个蒙特卡罗步骤，如下所示：
访问每个格点，
计算翻转自旋所需的能量：
dE如果其中是 [0,1] 中的随机数，则翻转自旋，否则离开它。
p
我从波动耗散定理中检索到的热容量：
我的代码重现了一些预期的行为：即在适当的温度下有一个关键的转变，但是：
能源有一个平滑的过渡。磁化是从 1 到 0 的急剧下降，而能量继续上升，并且没有显示出任何梯度变化。
热容量显示出一个峰值，但是当我增加晶格尺寸时，峰值并没有上升，而是变小了。
热容量很嘈杂：峰值经常翻倍，将其拟合到洛伦兹式是一场噩梦。更糟糕的是，多次运行模拟会得到相同的数据（如小误差条所示。
我的项目是一个模拟单个格子的 c++ 程序，以及一个分析输出的 python 脚本。到目前为止，脚本通过管道传输程序标准输出并将其分配给内部变量。它效率低下，但这是出于调试目的。
代码
完整的代码可以在这里找到：在此处输入链接描述。该代码是一个在单个格子上运行模拟的 c++ 程序和一个分析数据的 python 脚本。
重要的代码片段（即不包含不相关的内务处理）是
模拟.cpp
#include <iostream>
#include <math.h>
#include "include/simulation.h"
using namespace std;
/*
Advances the simulation a given number of steps, and updates/prints the statistics
into the given file pointer.
Defaults to stdout.
The number of time_steps is explcitly unsigned, so that linters/IDEs remind
the end user of the file that extra care needs to be taken, as well as to allow
advancing the simulation a larger number of times.
*/
void simulation::advance(unsigned int time_steps, FILE *output) {
unsigned int area = spin_lattice_.get_size() * spin_lattice_.get_size();
for (unsigned int i = 0; i < time_steps; i++) {
/*
To simulate an adiabatic Ising lattice, uncomment this
*/
// double temperature_delta = total_energy_/area - mean_energy_;
// if (abs(temperature_delta) < 1/area){
// cerr<<temperature_delta<<"! Reached equilibrium "<<endl;
// }
// temperature_ += temperature_delta;
if (time_ % print_interval_ == 0) {
total_magnetisation_ = spin_lattice_.total_magnetisation();
mean_magnetisation_ = total_magnetisation_ / area;
total_energy_ = compute_energy(spin_lattice_);
mean_energy_ = total_energy_ / area;
print_status(output);
}
advance();
}
}
/*
Advances the simulation a single step.
DOES NOT KEEP TRACK OF STATISTICS. Hence private.
*/
void simulation::advance() {
// #pragma omp parallel for collapse(2)
for (unsigned int row = 0; row < spin_lattice_.get_size(); row++) {
for (unsigned int col = 0; col < spin_lattice_.get_size(); col++) {
double dE = compute_dE(row, col);
double p = r_.random_uniform();
// float rnd = rand() / (RAND_MAX + 1.);
if (dE <0 || exp(-dE / temperature_) > p) {
spin_lattice_.flip(row, col);
}
}
}
time_++;
}
/*
Computes the total energy associated with spins in the spin_lattice_.
I originally used this function to test the code that tracked energy as the lattice
itself was modified, but that code turned out to be only marginally faster, and
not thread-safe. This is due to a race condition: when one thread uses a neighborhood
of a point, while another thread was computing the energy of one such point in
the neighborhood of (row, col).
*/
double simulation::compute_energy(lattice &other) {
double energy_sum = 0;
unsigned int max = other.get_size();
#pragma omp parallel for reduction(+ : energy_sum)
for (unsigned int i = 0; i < max; i++) {
for (unsigned int j = 0; j < max; j++) {
energy_sum += other.compute_point_energy(i, j);
}
}
return energy_sum/2;
}
void simulation::set_to_chequerboard(int step){
if (time_ !=0){
return;
}else{
for (unsigned int i=0; i< spin_lattice_.get_size(); ++i){
for (unsigned int j=0; j<spin_lattice_.get_size(); ++j){
if ((i/step)%2-(j/step)%2==0){
spin_lattice_.flip(i, j);
}
}
}
}
}
## 调查员.py
#! /usr/local/bin/python3
from os import mkdir
from os.path import exists
from subprocess import run, PIPE
from numpy import loadtxt, linspace, log, sqrt, append, exp, \
array, shape, inf, std, mean, argmax
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
class Simulation:
def __init__(self, magnetic_field=0.0, exchange_energy=1.0,
lattice_size=64, temperature=1.0, duration=50,
grain_size=0):
self.magnetic_field = magnetic_field
self.exchange_energy = exchange_energy
self.lattice_size = lattice_size
self.temperature = temperature
self.grain_size = grain_size
self.duration = duration
file_path = 'data/' + str(self) + '.csv'
if use_disk:
self.data = self.load_from_disk(file_path)
else:
self.data = self.run()
def times(self):
return self.data[:, 0]
def mean_magnetizations(self):
return self.data[:, 1]
def mean_energies(self):
return self.data[:, 2]
def load_from_disk(self, file_path):
if not exists(path=file_path):
self.run(file_path)
data = loadtxt(file_path)
data_duration = data[-1, 0] + self.duration/200 + 10
if self.duration > data_duration:
print('duration mismatch, ' + str(self.duration) + '!=' + str(
data_duration), end='. re-', flush=True)
self.run(file_path)
data = loadtxt(file_path)
return data
def run(self, file_path=None):
if file_path is None:
file_path = 'data/' + str(self) + '.csv'
# print('running simulation... ', self.duration, end=' ', flush=True)
command = ['./main', '-d', str(self.duration), '-t',
str(self.temperature), '-n', str(self.lattice_size),
'-j', str(self.exchange_energy), '-H',
str(self.magnetic_field), '-c', str(self.grain_size)]
if self.duration > 500:
command.append('-p')
command.append(str(int(self.duration/200)))
if use_disk:
command.append('-f')
command.append(file_path)
# print('Done! ')
run(command)
else:
r = run(command, stdout=PIPE)
raw_data = loadtxt(r.stdout.decode().split(sep='\n'),
delimiter='\t')
# print('Done! ')
return raw_data
def __str__(self):
return '(H=' + str(self.magnetic_field) + ')(J=' + str(
self.exchange_energy) + ')(T=' + str(
self.temperature) + ')(N=' + str(
self.lattice_size) + ')' + str(self.grain_size)
# ------------------------------------------------------------------------
def save_plot(title):
plt.legend(loc='best')
plt.suptitle(title)
file_name = title.replace(' ', '_')
if not exists('figures/'):
mkdir('figures/')
plt.savefig('figures/%s.png'%file_name)
plt.show()
# -------------------------------------------------------------------------
def smart_duration(temperature, multiplier=1.):
return int(((10**5)*base_duration*multiplier)/(
(temperature - t_c)**2 + breadth))
def theoretical_capacity(x, a, b, c, d):
return c + (b*x)/((x - a)**2 + d)
def investigate_heat_capacity(lattice_sizes=None, temps=None, **kwargs):
if temps is None:
temps = linspace(1.5, 3, 10)
if lattice_sizes is None:
lattice_sizes = [16,32, 34, 36]
fig, axes = plt.subplots(nrows=len(lattice_sizes), sharex='col')
crit_temps = []
for l, ax in zip(lattice_sizes, axes):
print('===',l)
crit_temps.append(fit_and_plot_capacity(ax, l, temps, **kwargs))
fig.set_size_inches(10.5, 10.5)
plt.xlabel('temperature / arb. u.')
save_plot('Heat capacity')
return crit_temps
def fit_and_plot_capacity(ax, l, temps, **kwargs):
"""
Plot the heat capacity of simulations at given temperature and lattice size.
Afterwards fit a lorentzian and plot.
Parameters:
-----------
ax : pyplot.axis
what to plot to.
l : int
lattice size
temps: numpy.array
Array of temperatures where to evaluate heat capacity.
Returns:
-------
popt[0]: float
most likely critical temperature.
"""
global use_disk
use_disk = False
simulations = [Simulation(lattice_size=l, temperature=t, duration=2)
for t in temps]
# sigmas = [stdev(s.mean_energies[:]) for s in simulations]
sigmas = array([multi_run(s, **kwargs) for s in simulations]).T[3]
meta_sigmas = array([multi_run(s, **kwargs) for s in simulations]).T[4]
# print(meta_sigmas)
Cs = [sigma**2/(temp**2)*10**3 for temp, sigma in zip(temps, sigmas)]
C_errs = Cs[:]*meta_sigmas[:]
try:
# popt, pcov = curve_fit(theoretical_capacity, temps*10, Cs,
# sigma=meta_sigmas, bounds=(
# [min(temps) - .2, 0, 0, 0],
# [max(temps) + .2, inf, inf, .7]))
pass
except RuntimeError:
print('I\'m too dumb to fit')
# popt = [temps[argmax(Cs)], (max(Cs) - min(Cs))/4, min(Cs),
# (max(temps) - min(temps))/4]
# ax.axvline(x=popt[0], ls='--', color='g')
# ax.plot(temps, theoretical_capacity(temps, *popt), 'g-',
# label='N = ' + str(l) + 'd = ' + str(popt[3]) + ' fit')
ax.errorbar(temps, Cs, fmt='b.', yerr=C_errs, label='N = ' + str(l))
ax.set_ylabel(r'C $\cdot 10^3$/ arb. u.')
ax.axvline(x=t_c, ls='-.', color='k')
ax.set_xticks(
append(linspace(min(temps), t_c, 6), linspace(t_c, max(temps), 6)))
ax.legend(loc='best')
# return popt[0]
def multi_run(sim, re_runs:int=1, take_last:int=300):
"""
Re run the Ising model simulation multiple times, and gather statistics.
Parameters:
----------
sim: simulation
re_runs: int
number of times to repeat simulation
take_last: int
How many of the final points to take statistics over.
Returns:
list:
"""
global use_disk
use_disk = False
sim.duration = smart_duration(sim.temperature)
print(sim.duration)
magnetizations = []
sigma_magnetizations = []
energies = []
sigma_energies = []
for i in range(re_runs):
sim.data = sim.run()
# Make each run take 100% longer, so that we
# can see if a system is still settling
last_magnetizations = sim.mean_magnetizations()[-take_last:]
magnetizations.append(abs(mean(last_magnetizations)))
sigma_magnetizations.append(std(last_magnetizations))
last_energies = sim.mean_energies()[-take_last:]
energies.append(mean(last_energies))
sigma_energies.append(std(last_energies))
return [mean(magnetizations), mean(sigma_magnetizations),
mean(energies), mean(sigma_energies),
std(sigma_energies)/mean(sigma_energies)]
# -------------------------------------------------------------------------
def finite_size_scale(N, t_inf, a, v):
return t_inf + a*(N**(-1/v))
def investigate_finite_size_scaling(critical_temperatures, lattice_sizes,
**kwargs):
if critical_temperatures is None:
critical_temperatures = investigate_heat_capacity(lattice_sizes,
**kwargs)
args, cov = curve_fit(finite_size_scale, lattice_sizes,
critical_temperatures)
plt.plot(lattice_sizes, critical_temperatures, 'b+', label='data')
plt.plot(lattice_sizes, finite_size_scale(lattice_sizes, *args), 'r-',
label='fit')
plt.ylabel('critical temperature / arb. u.')
plt.xlabel('Lattice size')
save_plot('Finite size scaling')
return args[0], sqrt(cov[0, 0])
# -------------------------------------------------------------------------
t_c = 2/log(1 + sqrt(2))
use_disk = False
breadth = 1
base_duration = 10
sizes = [16, 32]
# investigate_time_evolution()
# investigate_temperature_dependence(lattice_sizes=sizes, re_runs=4)
critical_temps = investigate_heat_capacity(lattice_sizes=sizes,
take_last=10000)
# temp_inf = investigate_finite_size_scaling(critical_temps, sizes)
# print(temp_inf)
# print((t_c - temp_inf[0])/ temp_inf[1], ' Standard errors away')
这是热容量，
和