参考文献
分子動力学法とモンテカルロ法
片岡洋右・著
発行 1994/05/10
準備
オンラインコンパイラを使用します。
プログラム
sample.py
import numpy as np
import math
# Constants
NSP = 2
EPS = 1e-2
SIGMA = 1e-0
TMASS = 1e-0
DT = 1e-3
NSTEP = 10000
# Initialize arrays
x = np.zeros(NSP)
y = np.zeros(NSP)
z = np.zeros(NSP)
xo = np.zeros(NSP)
yo = np.zeros(NSP)
zo = np.zeros(NSP)
xd = np.zeros(NSP)
yd = np.zeros(NSP)
zd = np.zeros(NSP)
vx = np.zeros(NSP)
vy = np.zeros(NSP)
vz = np.zeros(NSP)
fx = np.zeros(NSP)
fy = np.zeros(NSP)
fz = np.zeros(NSP)
# Initialize variables
ek = 0.0
ep = 0.0
h = 0.0
ho = 0.0
hsum = 0.0
h2sum = 0.0
eksum = 0.0
ek2sum = 0.0
sum = 0.0
def start():
print("MOLECULAR DYNAMICS SIMULATION BY VERLET METHOD")
print(f"NSP = {NSP}")
print(f"EPS = {EPS}, SIGMA = {SIGMA}")
print(f"TMASS = {TMASS}, DT = {DT}")
print(f"NSTEP = {NSTEP}")
def inconf():
global ek, ep, h, ho, hsum, h2sum, eksum, ek2sum, sum
xo[0] = -0.55
xo[1] = 0.55
for i in range(NSP):
yo[i] = 0.0
zo[i] = 0.0
xd[i] = 0.0
yd[i] = 0.0
zd[i] = 0.0
vx[i] = 0.0
vy[i] = 0.0
vz[i] = 0.0
ep = 0.0
for i in range(NSP - 1):
for j in range(i + 1, NSP):
rsq = (xo[j] - xo[i])**2 + (yo[j] - yo[i])**2 + (zo[j] - zo[i])**2
inv_rsqsq = 1.0 / rsq
inv_rsqsq3 = inv_rsqsq**3
inv_rsqsq6 = inv_rsqsq3**2
ep += 4.0 * EPS * (inv_rsqsq6**2 - inv_rsqsq6)
ek = 0.0
for i in range(NSP):
ek += 0.5 * TMASS * (vx[i]**2 + vy[i]**2 + vz[i]**2)
h = ek + ep
ho = h
hsum = ho
h2sum = ho**2
eksum = ek
ek2sum = ek**2
sum = 1.0
print(f"N = 0, HO = {ho:.6f}, EK = {ek:.6f}, EP = {ep:.6f}")
def force():
global ep
fx.fill(0.0)
fy.fill(0.0)
fz.fill(0.0)
ep = 0.0
for i in range(NSP):
for j in range(i + 1, NSP):
xx = (xd[j] - xd[i]) + (xo[j] - xo[i])
yy = (yd[j] - yd[i]) + (yo[j] - yo[i])
zz = (zd[j] - zd[i]) + (zo[j] - zo[i])
rsq = xx**2 + yy**2 + zz**2
if rsq != 0:
inv_rsqsq = 1.0 / rsq
inv_rsqsq3 = inv_rsqsq**3
inv_rsqsq6 = inv_rsqsq3**2
inv_rsqsq12 = inv_rsqsq6**2
ep_pair = 4.0 * EPS * (inv_rsqsq12 - inv_rsqsq6)
ep += ep_pair
dpdr = 4.0 * EPS * (-12.0 * inv_rsqsq12 + 6.0 * inv_rsqsq6) * inv_rsqsq
fx[i] += dpdr * xx
fy[i] += dpdr * yy
fz[i] += dpdr * zz
fx[j] -= dpdr * xx
fy[j] -= dpdr * yy
fz[j] -= dpdr * zz
def verlet(n):
global ek, h, hsum, h2sum, eksum, ek2sum, sum
old_fx = np.copy(fx)
old_fy = np.copy(fy)
old_fz = np.copy(fz)
if n == 0:
for i in range(NSP):
xd[i] = vx[i] * DT + 0.5 * fx[i] * DT**2 / TMASS
yd[i] = vy[i] * DT + 0.5 * fy[i] * DT**2 / TMASS
zd[i] = vz[i] * DT + 0.5 * fz[i] * DT**2 / TMASS
else:
for i in range(NSP):
xd[i] += fx[i] * DT**2 / TMASS
yd[i] += fy[i] * DT**2 / TMASS
zd[i] += fz[i] * DT**2 / TMASS
vx[i] += 0.5 * (old_fx[i] + fx[i]) * DT / TMASS
vy[i] += 0.5 * (old_fy[i] + fy[i]) * DT / TMASS
vz[i] += 0.5 * (old_fz[i] + fz[i]) * DT / TMASS
ek = 0.0
for i in range(NSP):
ek += 0.5 * TMASS * (vx[i]**2 + vy[i]**2 + vz[i]**2)
h = ek + ep
hsum += h
h2sum += h**2
eksum += ek
ek2sum += ek**2
sum += 1.0
def output(k, n):
if k == 0 and n % (NSTEP // 10) == 0:
print(f"N = {n}, H = {h:.6f}, EK = {ek:.6f}, EP = {ep:.6f}")
print(f"(H - HO) / HO = {(h - ho) / ho:.6f}, H - HO = {h - ho:.6f}")
if abs((h - ho) / ho) > 0.1:
print("CONSERVATION OF HAMILTONIAN IS POOR!!!!!!!!!!!!!!")
elif k == 1:
havg = hsum / sum
h2avg = h2sum / sum
ekavg = eksum / sum
ek2avg = ek2sum / sum
hrmd = math.sqrt(h2avg - havg**2)
ekrmd = math.sqrt(ek2avg - ekavg**2)
print(f"AVERAGE OF H = {havg:.6f}, RMD OF H = {hrmd:.6f}")
print(f"AVERAGE OF EK = {ekavg:.6f}, RMD OF EK = {ekrmd:.6f}")
print(f"NORMALIZATION CONSTANT = SUM = {sum:.6f}")
def main():
start()
n = 0
inconf()
force()
verlet(n)
output(0, n)
for n in range(1, NSTEP + 1):
force()
verlet(n)
output(0, n)
output(1, n)
main()
実行結果
console
MOLECULAR DYNAMICS SIMULATION BY VERLET METHOD
NSP = 2
EPS = 0.01, SIGMA = 1.0
TMASS = 1.0, DT = 0.001
NSTEP = 10000
N = 0, HO = -0.008684, EK = 0.000000, EP = -0.008684
N = 0, H = -0.008684, EK = 0.000000, EP = -0.008684
(H - HO) / HO = -0.000000, H - HO = 0.000000
N = 1000, H = -0.008051, EK = 0.000636, EP = -0.008687
(H - HO) / HO = -0.072907, H - HO = 0.000633
N = 2000, H = -0.006148, EK = 0.002542, EP = -0.008689
(H - HO) / HO = -0.292104, H - HO = 0.002537
CONSERVATION OF HAMILTONIAN IS POOR!!!!!!!!!!!!!!
N = 3000, H = -0.002975, EK = 0.005717, EP = -0.008692
(H - HO) / HO = -0.657428, H - HO = 0.005709
CONSERVATION OF HAMILTONIAN IS POOR!!!!!!!!!!!!!!
N = 4000, H = 0.001465, EK = 0.010160, EP = -0.008694
(H - HO) / HO = -1.168714, H - HO = 0.010149
CONSERVATION OF HAMILTONIAN IS POOR!!!!!!!!!!!!!!
N = 5000, H = 0.007171, EK = 0.015868, EP = -0.008697
(H - HO) / HO = -1.825797, H - HO = 0.015856
CONSERVATION OF HAMILTONIAN IS POOR!!!!!!!!!!!!!!
N = 6000, H = 0.014142, EK = 0.022842, EP = -0.008699
(H - HO) / HO = -2.628512, H - HO = 0.022827
CONSERVATION OF HAMILTONIAN IS POOR!!!!!!!!!!!!!!
N = 7000, H = 0.022377, EK = 0.031078, EP = -0.008702
(H - HO) / HO = -3.576691, H - HO = 0.031061
CONSERVATION OF HAMILTONIAN IS POOR!!!!!!!!!!!!!!
N = 8000, H = 0.031873, EK = 0.040577, EP = -0.008704
(H - HO) / HO = -4.670169, H - HO = 0.040557
CONSERVATION OF HAMILTONIAN IS POOR!!!!!!!!!!!!!!
N = 9000, H = 0.042629, EK = 0.051336, EP = -0.008707
(H - HO) / HO = -5.908777, H - HO = 0.051313
CONSERVATION OF HAMILTONIAN IS POOR!!!!!!!!!!!!!!
N = 10000, H = 0.054644, EK = 0.063354, EP = -0.008710
(H - HO) / HO = -7.292348, H - HO = 0.063328
CONSERVATION OF HAMILTONIAN IS POOR!!!!!!!!!!!!!!
AVERAGE OF H = 0.012440, RMD OF H = 0.018888
AVERAGE OF EK = 0.021137, RMD OF EK = 0.018895
NORMALIZATION CONSTANT = SUM = 10002.000000
[Execution complete with exit code 0]