基于Markov chains 计算一维谐振子基态能量

发表于

利用马可夫链计算一维谐振子的基态能量

计算作用量

$$ dS=\frac{(u_{n-1}+u_{n})^{2}}{2a\omega}+\frac{a\omega(u_{n-1}^{2}+u_{n}^{2})}{4} +\frac{(u_{n+1}+u_{n})^{2}}{2a\omega}+\frac{a\omega(u_{n+1}^{2}+u_{n}^{2})}{4} $$

在一系列旧点之上乱序选取，加上某个范围内的随机数，新点与旧点进行局部作用量比较，作用量如减小，则接受新点。如果作用量增大，则按照概率$e^{dS_{2}-dS{1}}$确定是否接受。不接受的点保留旧值。

1
2
3
4
5
double weight(double x0, double x1, double x2, double aomega){
	double S=0.0;
	S=(x1-x0)*(x1-x0)/(2.0*aomega) + aomega*(x0*x0+x1*x1)/4.0 + (x2-x1)*(x2-x1)/(2.0*aomega) + aomega*(x2*x2+x1*x1)/4.0;
	return exp(-1.0*S);
}

产生若干（很多）组态之后，对后面的组态求平方平均值即为基态能量。

1
2
3
4
5
import numpy as np
import matplotlib.pyplot as plt
import math
from ctypes import *
import numpy.ctypeslib as npct
1
2
libcd = npct.load_library("pf", ".")
array_1d_double = npct.ndpointer(dtype=np.float64, ndim=2, flags='CONTIGUOUS')
1
2
libcd.hhh.argtypes=[array_1d_double,c_int,c_int,c_double,c_double]
libcd.hhh.restype=c_int
1
lu=np.zeros([130000,500],dtype='float64')
1
libcd.hhh(lu,130000,500,0.1,0.2)
16
1
2
suu=np.zeros(130000)
sii=lu**2
1
2
for kk in range(0,130000):
    suu[kk]=np.mean(sii[kk,:])
1
plt.plot(suu[10000::600])#为防止关联性每600个取数据。
[<matplotlib.lines.Line2D at 0x7f8d80dda898>]

png

1
plt.plot(suu)
[<matplotlib.lines.Line2D at 0x7f8d405b7710>]

png

1
plt.hist(suu[100000::],40,histtype='bar',facecolor='yellowgreen',alpha=0.75) #方均值分布
(array([  59.,  204.,  421.,  739.,  653.,  868., 1014., 1371., 1563.,
        1890., 1927., 1802., 1762., 1831., 1505., 1209., 1271., 1342.,
        1111., 1041.,  907.,  809.,  658.,  616.,  520.,  426.,  294.,
         300.,  394.,  306.,  263.,  203.,  161.,  154.,  107.,   70.,
         109.,   86.,   27.,    7.]),
 array([0.30646848, 0.32209395, 0.33771941, 0.35334488, 0.36897034,
        0.3845958 , 0.40022127, 0.41584673, 0.4314722 , 0.44709766,
        0.46272313, 0.47834859, 0.49397406, 0.50959952, 0.52522499,
        0.54085045, 0.55647592, 0.57210138, 0.58772684, 0.60335231,
        0.61897777, 0.63460324, 0.6502287 , 0.66585417, 0.68147963,
        0.6971051 , 0.71273056, 0.72835603, 0.74398149, 0.75960695,
        0.77523242, 0.79085788, 0.80648335, 0.82210881, 0.83773428,
        0.85335974, 0.86898521, 0.88461067, 0.90023614, 0.9158616 ,
        0.93148707]),
 <a list of 40 Patch objects>)

png

1
np.mean(suu[10000::600])
0.49547505767135985
1
hy=lu[10000::600]
1
2
3
loo=np.zeros(500)
for po in range(0,500):
    loo[po]=np.mean(lu[10000,:]*lu[10000+po*1,:])
1
plt.plot((looo[5:50]),'rx')
[<matplotlib.lines.Line2D at 0x7f8d3f03e208>]

png

1
looo=np.log(loo)
1
loo
array([ 0.34280532,  0.13145158,  0.09545402,  0.05821956,  0.11446664,
       -0.00101133,  0.01600951,  0.01831628,  0.03161569,  0.03762436,
        0.04253305,  0.03503241, -0.01741445,  0.0206876 ,  0.0559742 ,
        0.00657246,  0.02660626,  0.04883725,  0.02803961,  0.04471572,
       -0.03246027,  0.02030995,  0.04995083,  0.08313026,  0.03374789,
        0.02868077, -0.04105371, -0.01369022, -0.02838291, -0.01871887,
        0.02223757,  0.04172896,  0.01848376,  0.06792485,  0.00939457,
       -0.0263474 , -0.00551463,  0.06894225, -0.01276136, -0.0205876 ])
1
np.var(suu[10000::600]) #f方差
0.010724637054091532

相关C++ 模块代码

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

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <iostream>
#include <fstream>
#include <random>
#include <vector>
#include <algorithm>

extern "C"
{
    int hhh(double * x,int hang,int lie, double aw, double dx );
}

double weight(double x0, double x1, double x2, double aomega){
	double S=0.0;
	S=(x1-x0)*(x1-x0)/(2.0*aomega) + aomega*(x0*x0+x1*x1)/4.0 + (x2-x1)*(x2-x1)/(2.0*aomega) + aomega*(x2*x2+x1*x1)/4.0;
	return exp(-1.0*S);
} 

int hhh(double * x,int hang,int lie, double aw, double dx ){
	int Nt,Nconf; 
	double aomega; // a(lattice spacing)*omega(angular frequency)
	double delta;  //update step size
	// char output_path[200]; // output file path
    aomega=aw;
    Nt=lie;
    Nconf=hang;
    aomega=aw;
    delta=dx;

// read in input
	// read_double("aomega", &aomega);
	// read_int("Nt", &Nt);
	// read_int("Nconf", &Nconf);
	// read_double("delta", &delta);
	// read_string("output_path", output_path);
	
	double* x1 = new double[Nt];
	double* x2 = new double[Nt];
	// double* x = new double[Nt*Nconf];
	std::ofstream out;

        std::default_random_engine seed;
       seed.seed(time(NULL));
        // seed.seed(1234);
    //    std::uniform_real_distribution<double> ran1(-0.1,0.1);
        std::uniform_real_distribution<double> ran2(-1.0,1.0);
        std::uniform_real_distribution<double> ran3(0,1.0);
	

// initialize the configuration to be a random number between [-0.1, 0.1]	
	for(int i=0; i<Nt; i++)
		// *(x1+i) = *(x2+i) = ran1(seed);
		*(x1+i) = *(x2+i) = 0.0;  // cold start

// update the configuration with metropolis algorithm
	std::vector<int> t_iter;
	for(int i=0; i<Nt; i++) t_iter.push_back(i);
	for(int n=0; n<Nconf; n++){ 
		std::random_shuffle (t_iter.begin(), t_iter.end());
		for(int t=0; t<Nt; t++){
			int i=t_iter[t];
			*(x2+i) = *(x1+i) + delta*ran2(seed);
			if(ran3(seed) < weight(*(x2+(i-1+Nt)%Nt),*(x2+i),*(x2+(i+1)%Nt),aomega)/weight(*(x1+(i-1+Nt)%Nt),*(x1+i),*(x1+(i+1)%Nt), aomega))
				*(x1+i)=*(x2+i);  // accept
			else
				*(x2+i)=*(x1+i);  // reject
		}
	
		for(int i=0; i<Nt; i++)
			*(x + n*Nt + i) = *(x2 + i);
	}

编译指令

1
2

g++ fp.cpp -shared -fPIC -o fp.so