Homec4science
Diffusion LibMultiScale (master) Browse

libmultiscale/test/single_scale/MD1D/test_1D_initial_temp/md1d.py
5228ea68cc05master

Last Change
Raw File

md1d.py
View Options

import numpy as np
import matplotlib.pyplot as plt
r0 = 1.23157195294235 # angstroms
# alpha = 0.0725672481645969
alpha = 0.
r0 = r0*(1+alpha)
print(r0)
NBATOMS = 99
mass = 39.95 # g/mol
timestep = 1 # femto-second
epsilon = 1.6567944e-21/6.94769e-21 # kcal/mol
sigma = 1.1
distance_standard = 1e-10
avogadro = 6.02214179e23
g2Kg = 1e-3
g_mol2Kg = g2Kg / avogadro
mass_standard = g_mol2Kg
boltzmann = 1.3806503e-23
cal2Joule = 4.184
kcal2Joule = 1e3 * cal2Joule
kcal_mol2Joule = kcal2Joule / avogadro
energy_standard = kcal_mol2Joule
time_standard = 1e-15
temperature_standard = 1.
force_standard = energy_standard / distance_standard
velocity_standard = distance_standard / time_standard
ft_m2v = (force_standard * time_standard / mass_standard) / velocity_standard
mvv2e = mass_standard * velocity_standard ** 2 / energy_standard
p0 = np.arange(NBATOMS)*r0
p = p0.copy()
v = np.zeros_like(p)
f = np.zeros_like(p)
a = np.zeros_like(p)
pe = np.zeros_like(p)
ke = np.zeros_like(p)
domainSize = p.shape[0] * r0
def computeForce():
pe[:] = 0
f[:] = 0
for i, x1 in enumerate(p):
for jj in [-2, 2, -1, 1]:
j = (i+jj) % p.shape[0]
x2 = p[j]
rij = x2 - x1
rij -= domainSize * round(rij / domainSize)
contribution = -24*epsilon/rij*(2*(sigma/rij)**12 - (sigma/rij)**6)
f[i] += contribution
epot = 4*epsilon*((sigma / rij)**12 - (sigma / rij)**6)
pe[i] += epot / 2
def linForce(r, r0):
f = 24*epsilon*sigma**6*(r0**6 - 2*sigma**6)/r0**13
f += 24*epsilon*sigma**6*(-7*r0**6 + 26*sigma**6)/r0**14*(r - r0)
return f
def quadEnergy(r, r0):
e = 4*epsilon*(-sigma**6/r0**6 + sigma**12/r0**12)
e += 24*epsilon*sigma**6*(r0**6 - 2*sigma**6)/r0**13*(r - r0)
e += 0.5*24*epsilon*sigma**6*(-7*r0**6 + 26*sigma**6)/r0**14*(r - r0)**2
return e
def computeLinearizedForce():
pe[:] = 0
f[:] = 0
for i, x1 in enumerate(p):
for jj in [-2, 2, -1, 1]:
j = (i+jj) % p.shape[0]
x2 = p[j]
rij = x2 - x1
rij -= domainSize * round(rij / domainSize)
# print(rij, jj*r0)
# assert np.fabs(rij - jj*r0) < 1e-1
contribution = linForce(rij, jj*r0)
f[i] += contribution
epot = quadEnergy(rij, jj*r0)
pe[i] += epot / 2
def computeEkin():
ke[:] = mvv2e * 0.5*mass*v**2
def computeTemperature():
computeEkin()
T = 2. / v.shape[0] / boltzmann * (ke.sum() * energy_standard)
return T
def computeStep(time_step):
computeEkin()
v[:] += 0.5 * time_step * a
p[:] += time_step * v
computeForce()
# computeLinearizedForce()
a[:] = ft_m2v * f / mass
v[:] += 0.5 * time_step * a
def plot_disp(u):
plt.plot(u, label='disp')
plt.show()
epot = []
ekin = []
# computeForce()
computeLinearizedForce()
ref_epot = pe.sum()
print(ref_epot)
def rescaleV(temp):
print(temp)
T = computeTemperature()
v[:] *= np.sqrt(temp/T)
computeEkin()
T = computeTemperature()
print(T)
plot_disp(v)
def initT(temp):
# v[:] = np.random.normal(loc=0, size=v.shape[0])
v[:] = np.random.random(size=v.shape[0])
v[:] -= np.average(v)
rescaleV(temp)
initT(100)
# p[15] += 1e-4
T = []
disp = []
for i in range(1000):
computeStep(timestep)
t = computeTemperature()
if (i % 100 == 0):
print(i, t)
epot.append(pe.sum())
ekin.append(ke.sum())
T.append(t)
if (i % 10 == 0):
# plot_disp(p - p0)
disp.append(p - p0)
disp = np.array(disp)
epot = np.array(epot) # - ref_epot
ekin = np.array(ekin)
etot = epot+ekin
plt.plot(epot, label='epot')
plt.plot(ekin, label='ekin')
plt.plot(etot, label='etot')
plt.legend()
plt.show()
plt.plot(T, label='temperature')
plt.legend()
plt.show()
plt.imshow(disp)
plt.show()
plot_disp(disp[-1, :])
print(np.average(T[100:]))
ekin = ekin[100:]
epot = epot[100:]
etot = etot[100:]
print(etot.min(), etot.max(), (etot.max() - etot.min())/etot.max())
print(ekin.min(), ekin.max(), (ekin.max() - ekin.min())/ekin.max())
print(epot.min(), epot.max(), (epot.max() - epot.min())/epot.max())
c4science · Help