计算作业2:理想气体中的径向分布

理想气体中的径向分布

首先，我们必须生成理想气体的“configuration”,我们将在二维平面做。这相对更加容易，因为它仅表示每个粒子的x和y在0和L之间随机。

1
2


import numpy as np
import matplotlib.pyplot as plt

1.设置生成粒子

 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


# 返回生成的设置
def generate_configuration(N_particles,box_L):
    '''
    :param N_particles:int,粒子数
    :param box_L: float,盒子长度
    :return: numpy.array
    '''
    # 使用np.random.randon进行生成，生成的array，N_particles作为行，x,y作为两列，xy的位置位于0到box_L

    configuration = box_L * np.random.random(size=(N_particles,2))

    return configuration

# 绘画
def plot_configuration(configuration,box_L):
    plt.figure(figsize=(4,4))
    plt.axis('equal')
    plt.scatter(configuration[:,0],configuration[:,1])
    plt.xlim(0,box_L)
    plt.ylim(0,box_L)


#测试是否正常
test_box_L = 20
test_num_particles = 100
configuration1 = generate_configuration(test_num_particles,test_box_L)
plot_configuration(configuration1,test_box_L)

png

2.生成距离列表

 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


# 使用一个函数来获取两两粒子之间的距离
def compute_distances(configuration):
    distance_list = []
    num_particles = configuration.shape[0]
    for i in range(num_particles):
        for j in range(num_particles):
            #排除粒子自身距离的计算
            if i == j :continue
            posi = configuration[i]
            posj = configuration[j]
            #我们计算pos1和pos2之间的欧几里得距离 设置变量为dist

            # dr 为 向量 dx,dy
            dr = posj -posi
            #dr2 为向量 (dx*dx,dy*dy)
            dr2 = dr*dr
            # dist = sqrt(dx^2 + dy^2)
            dist = np.sqrt(dr2.sum())

            distance_list.append(dist)
    return np.array(distance_list)

# 返回距离列表的长度应为n*（n-1），让我们检查一下
distance_list_1 = compute_distances(configuration1)
print("Does the number of distances equal the number expected?")
print(len(distance_list_1),"=", test_num_particles*(test_num_particles-1) )

Does the number of distances equal the number expected?
9900 = 9900

3.用柱状图展示距离

现在我们用直方图来展示距离列表，这意味着会计算r到r+dr之间的频率，可以使用numpy的histogram函数进行计算。我们将会进行计算0-box_L/2之间的统计，具体原因会在第4部分解释。 numpy histogram的详细文档可以参考: 官方文档

1
2
3


# 首先让我们检查粒子之间的距离
print("Min distance:",distance_list_1.min())
print("Max distance:",distance_list_1.max())

Min distance: 0.055340547910658046
Max distance: 26.320208630702233

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12


def histogram_distances_default(distance_list):
    hist, bin_edges = np.histogram( distance_list )
    return hist, bin_edges

def plot_histogram(hist,bin_edges):
    bin_centers = (bin_edges[:-1]+bin_edges[1:])/2.0
    plt.plot(bin_centers,hist,marker='o')
    plt.ylabel("N(r)")
    plt.xlabel("$r$")

dist_hist_1, bin_edges_1 = histogram_distances_default(distance_list_1)
plot_histogram(dist_hist_1,bin_edges_1)

png

从图上可以看出，r到r+dr距离的粒子先是增加，随后减少，这就是边际效应，因为会达到边缘。两个粒子的最大距离是L*sqrt(2),Why?

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15


def histogram_distances(distance_list, max_dist, bin_size):
    # 我们修改一下函数，能够控制多少bins被计算（柱状图柱子数量）
    bins = np.arange(0, max_dist+bin_size, bin_size)
    hist, bin_edges = np.histogram( distance_list, bins=bins )
    return hist, bin_edges


#我们选择一个箱子数让其看起来更加平滑
#设置最大的距离为L*sqrt(2)/2

test_bin_size = 0.5
max_dist = test_box_L/2.0*np.sqrt(2)
dist_hist_1, bin_edges_1 = histogram_distances(distance_list_1,max_dist=max_dist, \
                                               bin_size=test_bin_size)
plot_histogram(dist_hist_1,bin_edges_1)

png

4. 周期性条件

在模拟中，|我们始终具有有限系统大小的限制。为了摆脱一些边界效应，我们可以使用周期性的边界条件。这意味着我们计算距离，好像边界像甜甜圈一样缠绕。我们使用的距离计算称为最小图像约定。这意味着两个粒子之间的距离是系统无限个副本中两个粒子之间的最小距离。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15


def plot_configuration_pbc(configuration,box_L):
    plt.figure(figsize=(4,4))
    plt.axis('equal')
    for shift_x in range(-1,2):
        for shift_y in range(-1,2):
            if shift_x == 0 and shift_y == 0:
                alpha = 1
            else:
                alpha = .3
            plt.scatter(shift_x*box_L+configuration[:,0],shift_y*box_L+configuration[:,1],alpha=alpha)

    plt.xlim(-box_L,2*box_L)
    plt.ylim(-box_L,2*box_L)

plot_configuration_pbc(configuration1,test_box_L)

png

我们重新写一个最小图像距离计算器

你可以和之前的对比看有什么不同之处

 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


def compute_distances_minimum_image( configuration, box_size ):
    distance_list = []
    num_particles = configuration.shape[0]
    for i in range(num_particles):
        for j in range(num_particles):
            if i == j: continue
            posi = configuration[i]
            posj = configuration[j]

            dr = posj-posi
            #不同之处
            #最小图像计算的 dr - 你知道这么计算的原因吗?
            dr = dr - box_size*np.floor(dr/box_size+0.5)

            dr2 = dr*dr
            dist = np.sqrt( dr2.sum() )
            distance_list.append(dist)

    return np.array(distance_list)

distance_list_1_pbc = compute_distances_minimum_image(configuration1, box_size = test_box_L)

#现在没有两个粒子是超过 L*sqrt(2)/2 这个距离的，你知道为什么吗?
print("Min distance:",distance_list_1_pbc.min())
print("Max distance:",distance_list_1_pbc.max())

test_bin_size = .25
#将图像界面设置最大距离为 L*sqrt(2)/2
max_dist = test_box_L/2.0*np.sqrt(2)
dist_hist_1_pbc, bin_edges_1_pbc = histogram_distances(distance_list_1_pbc,max_dist=max_dist, bin_size=test_bin_size)
plot_histogram(dist_hist_1_pbc,bin_edges_1_pbc)
plt.axvline(test_box_L/2.,color='black',linestyle='--')

Min distance: 0.055340547910658046
Max distance: 14.0151727079517





<matplotlib.lines.Line2D at 0x1fa54ed3760>

png

您应该看到在L/2的该图中发生一些有趣的事情。尝试想象上面显示的系统的3x3副本，然后查看哪些粒子在半径为r的圆内。您在上面的图中看到半径为r的圆和半径为r的正方形之间的面积差有关。因此我们只需要查看L/2的内容？？？？ Why？？？随后我们重新绘制到距离到 L/2并且重新选择一个bins，使得图像更加光滑

5.计算径向分布函数

径向分布函数的定义为该直方图除以该点半径r内部点的数量。在二维中，为2 pi r dr * density 在之前我们计算了每个粒子的距离，这意味着基本上为每个粒子计算了一个直方图并且将它们组合了起来，因此我们还要除以N 在get_gofr函数中，我们计算bin_centers数组和数的密度。最终计算径向分布函数。储存在gofr的数组中，并将其绘画出来。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17


def plot_rdf(gofr,bin_centers):
    plt.plot(bin_centers,gofr,marker='o')
    plt.ylabel("g(r)")
    plt.xlabel("$r$")

def get_gofr(hist,bin_edges,num_particles, box_size):
    rho = num_particles/box_size/box_size
    bin_centers = (bin_edges[1:]+bin_edges[:-1])/2.0
    dr = bin_edges[1]-bin_edges[0]
    denominator = 2.*np.pi*bin_centers*dr*rho*( num_particles )
    gofr = hist/denominator

    return gofr, bin_centers

gofr, bin_centers = get_gofr( dist_hist_1_pbc, bin_edges_1_pbc, test_num_particles, test_box_L )

plot_rdf(gofr, bin_centers)

png

平均径向分布

当我们进行模拟的时候，我们会生成许多设置，求径向分布函数的平均值。现在让我们对我们理想的气体做这件事。在下面的代码中我们生成N个配置，计算他们的g(r)使用上面的函数，增加num_configurations知道g(r)的误差每个点小于5%

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17


num_configurations = 50
gofr_list = []
for i in range(num_configurations):
    #生成配置
    configuration_i = generate_configuration(test_num_particles,test_box_L)
    #计算pbc距离
    distance_list_i_pbc = compute_distances_minimum_image(configuration_i, box_size = test_box_L)
    #距离直方图
    dist_hist_i_pbc, bin_edges_i_pbc = histogram_distances(distance_list_i_pbc,max_dist=test_box_L/2., bin_size=test_bin_size)
    # get_gofr
    gofr, bin_centers = get_gofr( dist_hist_i_pbc, bin_edges_i_pbc, test_num_particles, test_box_L )
    gofr_list.append(gofr)

avg_gofr = np.mean(gofr_list, axis=0)
plot_rdf(avg_gofr, bin_centers)
plt.axhline(1.0,linestyle='--',color='black')
plt.ylim(0.8,1.2)

(0.8, 1.2)

png