-
Notifications
You must be signed in to change notification settings - Fork 2
/
force_field_hexagon.py
150 lines (120 loc) · 3.1 KB
/
force_field_hexagon.py
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
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
#coding=utf-8
import math
import random as random1
import sys
import numpy as np
def to_po(num):
'''将随机选取的序数变为坐标'''
po=[]
po_cen=[]
for i in range(3):
for j in range(3):
po.append([num/n+n*i,num%n+n*j])
po_cen.append([num/n+n*1,num%n+n*1])
return(po,po_cen)
def po_sele(seq,sel):
'''随机选取原子位置'''
sele_num=[]
sele_po=[]
sel_po_cen=[]
while (len(sele_num) < sel):
sele_num.append(random1.choice(seq))
sele_num = list(set(sele_num))
for i in sele_num:
(po,po_cen)=to_po(i)
sele_po.extend(po)
sel_po_cen.extend(po_cen)
return(sele_po,sele_num,sel_po_cen)
def hex_en(po,sele_po):
'''用拟合关系计算一个原子的能量'''
x,y=po[0],po[1]
hexagon=[[x,y-1],[x+1,y-1],[x+1,y],[x,y+1],[x-1,y+1],[x-1,y]]
hex_po=[]
dis_hex_po=[]
for i in hexagon :
if i in sele_po :
hex_po.append(i)
len_hex_po = len(hex_po)
if len_hex_po ==0 :
return force_field[0]
elif len_hex_po ==1 :
return force_field[1]
elif len_hex_po ==5 :
return force_field[11]
elif len_hex_po ==6 :
return force_field[12]
else :
for i in hex_po :
for j in hex_po :
dis=np.sum((a[i[0],i[1],]-a[j[0],j[1],])**2)
dis_hex_po.append(round(dis,0))
diff_13 = dis_hex_po.count(1)-dis_hex_po.count(3)
dis_hex_po = sorted(set(dis_hex_po))
if len_hex_po == 2:
if dis_hex_po == [0.0, 1.0] :#临
return force_field[2]
elif dis_hex_po == [0.0, 3.0] :#间
return force_field[3]
elif dis_hex_po == [0.0, 4.0] :#对
return force_field[4]
if len_hex_po == 3:
if dis_hex_po == [0.0, 1.0, 3.0] : #连
return force_field[5]
elif dis_hex_po == [0.0, 1.0, 3.0, 4.0]: #偏
return force_field[6]
elif dis_hex_po == [0.0, 3.0] : #均
return force_field[7]
if len_hex_po == 4:
if diff_13 == 2 :
return force_field[8]
elif diff_13 == -2:
return force_field[9]
elif diff_13 == 0:
return force_field[10]
def re_po_sele(sele_num):
'''改变一个原子位置 计算能量'''
b=sele_num[0]
while (b in sele_num):
b=random1.choice(range(n**2))
sele_num.pop(random1.choice(range(sel)))
sele_num.append(b)
sel_po,sel_po_cen=[],[]
for i in sele_num:
(po,po_cen)=to_po(i)
sele_po.extend(po)
sel_po_cen.extend(po_cen)
#return(sele_po,sele_num,sel_po_cen)
energy=0
for i in sele_po_cen :
energy=energy+hex_en(i,sele_po)
return(energy,sele_num)
def b_z(num,sele_num):
(energy,sele_num)=re_po_sele(sele_num)
(energy1,sele_num1)=re_po_sele(sele_num)
prob=np.exp((energy-energy1)/kt)
txt=open(str(sel)+"_in_"+str(n),"a")
i=0
while (i<num):
#print prob
if prob > random1.uniform(0,1) :
txt.write(str(energy1)+" "+str(prob)+"\n")
energy,sele_num=energy1,sele_num1
else:
txt.write(str(energy)+" "+str(prob)+"\n")
energy1,sele_num1=re_po_sele(sele_num)
prob=np.exp((energy-energy1)/kt)
i=i+1
n=3
sel=5
force_field=range(13)
a = np.linspace(0,0,2*(n*3)**2).reshape(n*3,n*3,2)
a.dtype=np.float64
#print a.shape
for i in range(n*3):
for j in range(n*3):
a[i,j,0] = i+0.50*j
a[i,j,1] = j*math.sqrt(3.0)*0.5
(sele_po,sele_num,sele_po_cen)=po_sele(range(n**2),sel)
print hex_en([3,3],sele_po)
print sele_num
b_z(1000,sele_num)