Animation が動作しない

前提・実現したいこと

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_yiy

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/ツールのバージョンなど）

ここにより詳細な情報を記載してください。