-
Notifications
You must be signed in to change notification settings - Fork 0
/
solution.visual.02.py
167 lines (119 loc) · 5.47 KB
/
solution.visual.02.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
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
# Граничные условия (отражающая стенка)
# IMPORTS
import numpy as np
from vpython import *
# CONSTANTS
SIGMA = 0.272
BOILING_POINT = 27.1 # https://en.wikipedia.org/wiki/Neon
MELTING_POINT = 24.55
TEMP_EPS = 47.0
TEMP = BOILING_POINT / TEMP_EPS # температура релаксации
TEMP_0 = MELTING_POINT / TEMP_EPS
TEMP_1 = 2 * TEMP
CELL_SIZE = 1.2 * SIGMA
RANDOM_AREA = 0.8 * SIGMA
CELL_IDENT = (CELL_SIZE - RANDOM_AREA) / 2
GRID_SIZE = 10
COUB_SIZE = CELL_SIZE * GRID_SIZE
N_ATOM = GRID_SIZE ** 3
ATOM_RADIUS = CELL_SIZE / 10
CAMERA_SIZE = 500
CAMERA_POS = CELL_SIZE * GRID_SIZE / 2
FULL_TIME = 3e-5
INTEGRATION_STEP = 3e-8
ITERATIONS = int(FULL_TIME / INTEGRATION_STEP)
# FUNCTIONS
def calcPosByCell(cellPos):
""" Рассчитывает координаты атома по положению ячейки """
initCoord = CELL_SIZE * cellPos
identArray = CELL_IDENT * np.ones(3)
randomArray = np.random.uniform(0, RANDOM_AREA, 3)
return initCoord + identArray + randomArray
def initPositions():
""" Строит первоначальное расположение атомов """
initPos = np.empty(N_ATOM * 3).reshape(N_ATOM, 3) # Выделяем память
t = 0
for x in range(GRID_SIZE):
for y in range(GRID_SIZE):
for z in range(GRID_SIZE):
initPos[t] = calcPosByCell(np.array([x, y, z]))
t += 1
return initPos
def initAcceleration():
""" Строит нулевой массив ускорений """
return np.zeros(N_ATOM*3).reshape(N_ATOM, 3)
def initSpeed():
""" Строит случайные начальные скорости """
v_ = np.random.uniform(-1, 1, (N_ATOM, 3))
# нормировка амплитуд скоростей, чтобы кинетическая энергия соотвествовала начальной температуре
velsq = np.sum(v_ ** 2)
aheat = 3.0 * N_ATOM * INTEGRATION_STEP ** 2 * TEMP
factor = np.sqrt(aheat / velsq)
return v_ * factor
def initForces():
""" Строит нулевой массив векторов сил """
return np.zeros(N_ATOM*3).reshape(N_ATOM, 3)
def evaluateForces(coords, forces, energy = 0):
""" Вычисление сил действующих на частицы | вычисление энергии"""
for i in range(N_ATOM - 1):
# Вычисляется разность между радиус векторами частиц
ijCoord = coords[i] - coords[i+1:]
# Вычисляем вспомог расстояния между атомами
dist2 = np.sum(ijCoord * ijCoord, axis=1)
dist2Inv = 1.0 / dist2
dist6Inv = dist2Inv * dist2Inv * dist2Inv
# Вычисляем межатомные силы
ff = 48.0 * dist2Inv * dist6Inv * (dist6Inv - 0.5)
force_ = np.einsum('ki,k->ki', ijCoord, ff)
forces[i] += np.sum(force_, axis=0)
forces[i+1:] -= force_
# Потенциальная энергия Леннарда-Джонса
enr = 4.0 * dist6Inv * (dist6Inv - 1.0)
energy += np.sum(enr)
return energy
def werleScheme(acceleration, speed, coords, forces):
""" Решение уравнений движения по 3й схеме Верле """
# acceleration = forces * INTEGRATION_STEP * INTEGRATION_STEP / 2
# speed += 2.0 * acceleration
speed += (acceleration + forces) / 2 * INTEGRATION_STEP
acceleration = forces
forces[:, :] = 0
coords += speed * INTEGRATION_STEP + acceleration * INTEGRATION_STEP * INTEGRATION_STEP / 2
# VPYTHON FUNCTIONS
def createAtomsByPos(pos):
""" Создание сфер-атомов """
atoms = []
for i in range(len(pos)):
P = vector(pos[i][0], pos[i][1], pos[i][2])
atoms.append(sphere(pos=P, radius=ATOM_RADIUS))
return atoms
def changePosAtoms(atoms, positions):
""" Изменение позиций сфер-атомов """
for i in range(len(atoms)):
P = positions[i]
atoms[i].pos = vector(P[0], P[1], P[2])
# PREPARATION
scene = canvas(width=CAMERA_SIZE, height=CAMERA_SIZE)
scene.camera.pos = vector(CAMERA_POS, CAMERA_POS, 0)
scene.autoscale = False
graph(width = 600, height = 200, title = 'Total energy of the system')
curve = gcurve(color=color.orange)
coords = initPositions() # массив радиус-векторов частиц
acceleration = initAcceleration() # массив ускорений
speed = initSpeed() # массив скоростей
forces = initForces() # массив векторов сил, действующих на каждую частицу
atoms = createAtomsByPos(coords)
# LIFECYCLE
for iter in range(ITERATIONS):
""" Вычисление сил действующих на частицы + потенциальная энергия """
energy = evaluateForces(coords, forces)
""" Решение уравнений движения / 3 схема Верле"""
werleScheme(acceleration, speed, coords, forces)
# Отталкивание c изменением температуры
for i in range(N_ATOM - 1):
for j in range(3):
if (coords[i][j] <= 0) or (coords[i][j] >= COUB_SIZE):
speed[i] *= -np.sqrt(TEMP_0 / TEMP)
changePosAtoms(atoms, coords)
energy += np.sum(speed ** 2)
curve.plot(pos=(iter, energy))