質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

新規登録して質問してみよう
ただいま回答率
85.35%
Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

Q&A

0回答

419閲覧

Animation が動作しない

HIROKI123

総合スコア0

Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

0グッド

0クリップ

投稿2020/10/27 06:27

前提・実現したいこと

Pythonで分子シミュレーションをしたく、各サイトにあるものをいくつか回していました。
そこで、プロットを動画にしたいと思いAnimationを用いて描画したいと思っています。
https://qiita.com/sci_Haru/items/2c8874c26bedab04088a#%E5%86%85%E5%AE%B9

発生している問題・エラーメッセージ

Figureそのものは表示されるのですが、プロットや軸が全く表示されません。

該当のソースコード

import numpy as np
import random
import matplotlib.pyplot as plt
import matplotlib.animation as animation

fig = plt.figure()
ims = []

lat_x = 0.691 # 格子定数: x
lat_y = 0.691
Lx = 5 # 格子点: x
Ly = 5
N_atom = Lx*Ly
m=1 # 質量

delta_t = 0.002

T_const = 1 # 無次元量 #T_const = T_const /119 # convert to natural unit

Max_step =100

x = np.zeros([N_atom],float)
y = np.zeros([N_atom],float)

vx = np.zeros([N_atom],float)
vy = np.zeros([N_atom],float)

fx = np.zeros([N_atom,2],float)# Q: What is this "2"?
fy = np.zeros([N_atom,2],float)

T_lis = []# step毎の物理量の格納
P_lis = []
KE_lis =[]
PE_lis =[]
Eot_lis =[]
Tav_lis = []
Pav_lis = []
KEav_lis =[]
PEav_lis =[]
Eotav_lis =[]

density = N_atom/(LxLylat_x*lat_y) # 数密度: 2次元

print ("number density = ", density)

def Maxwell_velocity2D(TT): # マクスウェルの速度分布から: ボックス-ミュラー法: self
for i in range(N_atom):
R1=random.random()
R2=random.random()
R3=random.random()
R4=random.random()
vx[i] = (np.sqrt(-2* (TT/m)np.log(R1)))np.cos(2np.piR2)
vy[i] = (np.sqrt(-2* (TT/m)np.log(R3)))np.cos(2np.piR4)

def initialposMaxwellvel():
i=-1
for ix in range(Lx):
for iy in range(Ly):
i += 1
x[i] = lat_xix
y[i] = lat_y
iy

Maxwell_velocity2D(T_const)

def velocity_scaling(KE): #速度スケーリング法による温度制御

fac_T=np.sqrt(N_atom*T_const/(KE)) vx[:] = fac_T*vx[:] vy[:] = fac_T*vy[:]

def sign(a,b):
if b >= 0:
return abs(a)
else:
return -abs(a)

def Forces(t,w,PE,PEorW):
rcut = 4#カットオフ
PE = 0

r2cut = rcut**2 fx[:,t] = 0 fy[:,t] = 0 for i in range(N_atom-1): for j in range(i+1, N_atom): dx = x[i] - x[j] dy = y[i] - y[j] if (abs(dx) > 0.5*lat_x*Lx): # 近接 イメージ相互作用 dx = dx - sign(lat_x*Lx,dx) if (abs(dy) > 0.5*lat_y*Ly): dy = dy-sign(lat_y*Ly,dy) r2 = dx**2+dy**2 if (r2 < r2cut): if r2 == 0 : r2 = 0.0001 invr2 = 1/r2 wij = 48*(invr2**3-0.5)*(invr2**3) fijx = wij*invr2*dx fijy = wij*invr2*dy fx[i,t] += fijx fy[i,t] += fijy fx[j,t] -= fijx fy[j,t] -= fijy PE += 4*(invr2**3)*((invr2**3)-1) w += wij if (PEorW == 1) : return PE else : return w

def time_evolution():# メイン: 時間発展

avT = 0 avP = 0 Pavg = 0 avKE=0 avPE=0 t1 = 0 PE = 0.0 KE = 0 w = 0 initialposMaxwellvel() KE = (np.sum(vx**2)+np.sum(vy**2))/2 velocity_scaling(KE) #速度スケーリングによる温度一定化 KE = (np.sum(vx**2)+np.sum(vy**2))/2 T = KE/N_atom P = density*(2*KE+1.5*w)/(3*N_atom) # 圧力の計算: for 2D PE = Forces(t1, w, PE, 1) time = 1 n_step=0 while n_step < Max_step: if n_step % 50 ==0: print ("step=",n_step) print("T=",T," KE/N = ",KE/N_atom," PE/N = ",PE/N_atom," Etot/N = ",(KE+PE)/N_atom, " P=",P) if n_step % 5 ==0: im=plt.scatter(x, y, s=500, c="pink", alpha=0.9, linewidths="1.5",edgecolors="black") ims.append([im]) n_step += 1 for i in range(N_atom): # PE = Forces(t1, w, PE, 1) x[i] +=delta_t*(vx[i]+delta_t*fx[i,t1]/2) # 速度ベルレー法による更新 y[i] +=delta_t*(vy[i]+delta_t*fy[i,t1]/2) if x[i] <= 0: x[i] += lat_x*Lx if x[i] >= lat_x*Lx: x[i] -= lat_x*Lx if y[i] <= 0: y[i] += lat_y*Ly if y[i] >= lat_y*Ly: y[i] -= lat_y*Ly PE = 0. t2=1 PE = Forces(t2, w, PE, 1) KE = 0. w = 0. vx[:] += delta_t*(fx[:,t1]+fx[:,t2])/2 # 速度ベルレー法による更新 vy[:] += delta_t*(fy[:,t1]+fy[:,t2])/2 KE =(np.sum(vx**2)+np.sum(vy**2))/2 w = Forces(t2, w, PE, 2) velocity_scaling(KE) KE = (np.sum(vx**2)+np.sum(vy**2))/2 T = KE/N_atom P = density*(2*KE+1.5*w)/(3*N_atom) # for 2D avT += T avP += P avKE += KE avPE += PE time +=1 t = time if (t == 0): t = 1 Pavg =avP/t eKavg = avKE/t ePavg = avPE/t Tavg = avT/t T_lis.append(T) # P_lis.append(P) # P KE_lis.append(KE/N_atom) PE_lis.append(PE/N_atom) dumEtot = KE+PE Eot_lis.append( dumEtot /N_atom) Tav_lis.append(Tavg )# Pav_lis.append(Pavg) KEav_lis.append(eKavg/N_atom) PEav_lis.append(ePavg/N_atom) Eotav_lis.append((eKavg+ePavg)/N_atom)

time_evolution()

plt.axis([0,lat_xLx,0,lat_yLy])
ani = animation.ArtistAnimation(fig, ims, interval =1)
plt.show()
ani.save("output.gif", writer="imagemagick")

### 試したこと Code内のlistに問題があるのかいくつか検討をしました。 また、他のCodeでアニメーションを試したところ動作したので、環境的には問題ないと思ています。 https://qiita.com/yubais/items/c95ba9ff1b23dd33fde2 ### 補足情報(FW/ツールのバージョンなど) ここにより詳細な情報を記載してください。

気になる質問をクリップする

クリップした質問は、後からいつでもMYページで確認できます。

またクリップした質問に回答があった際、通知やメールを受け取ることができます。

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

jeanbiego

2020/10/27 06:32

このままだとインデントわかりませんので、コードは```と```の間に書くようにしてください。(しようとした形跡はありますが…)
guest

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

まだ回答がついていません

会員登録して回答してみよう

アカウントをお持ちの方は

15分調べてもわからないことは
teratailで質問しよう!

ただいまの回答率
85.35%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

関連した質問