Python分子动力学

1. 前言
- 1.1 分子热运动
- 1.2 理想气体模型
- - 1.2.1 弹性碰撞
  - 1.2.2 二维分子体系的弹性碰撞
- 1.3 理想气体分子动力学模拟
- 1.4 Lennard-Jones potential
- 1.5 Periodic boundary conditions
- 1.6 Reduced units
- 1.7 Calculating the forces
- 1.8 Temperature
- 1.9 Molecular dynamics program
- - 1.9.1 Initialising the particles
  - 1.9.2 Integrating the equations of motion

1. 前言

分子动力学是一套分子模拟方法，该方法主要是依靠计算机来模拟分子、原子体系的运动，是一种多体模拟方法。通过对分子、原子在一定时间内运动状态的模拟，从而以动态观点考察系统随时间演化的行为。通常，分子、原子的轨迹是通过数值求解牛顿运动方程得到，势能（或其对笛卡尔坐标的一阶偏导数，即力）通常可以由分子间相互作用势能函数、分子力学力场、全始计算给出。对于考虑分子本身的量子效应的体系，往往采用波包近似处理或采用量子力学的费恩曼路径积分表述方式[1]处理。分子动力学也常常被采用作为研究复杂体系热力学性质的采样方法。在分子体系的不同状态构成的系综中抽取样本，从而计算体系的构型积分，并以构型积分的结果为基础进一步计算体系的热力学量和其他宏观性质。分子动力学最早在20世纪50年代由物理学家提出，如今广泛应用于物理、化学、生物体系的理论研究中。

1.1 分子热运动

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-ikwlSeKR-1654083484608)(./Images/Translational_motion.gif)]

上图为一个分子系统热运动的模拟动画。其中一些小球被涂成了红色方便追踪。在这个模拟中，平均而言，两个分子碰撞前后动能保持不变。同时，我们忽视了黑体辐射，即系统的总能量保持不变。

1.2 理想气体模型

1.2.1 弹性碰撞

液体或者气体中分子之间发生的碰撞并不是完美的弹性碰撞。因为在碰撞的过程中，分子的平动会和其他运动自由度之间发生能量交换。但是平均而言，分子间的碰撞可以被简化为弹性碰撞，即分子间只发生平动动能交换。

我们考虑如下一个简化系统作为分子热运动的近似。

分子抽象为一个具有质量m和直径r的圆球。
分子间的相互作用力暂时被忽略。
当两个分子接触时，发生完全弹性碰撞。

1.2.2 二维分子体系的弹性碰撞

为了方便理解和作图，我们这里对二维情况下的弹性碰撞进行计算。对于两个无旋转的刚体而言，发生弹性碰撞时遵从动量守恒，能量守恒和角动量守恒。如果二者之间没有摩擦力，则二者速度在碰撞点切线方向上保持不变，碰撞只改变速度在垂直切线上的分量。

碰撞之后二者速度的表达式为：

上式中括号表示两个矢量的内积或者点乘。

1.3 理想气体分子动力学模拟

我们将使用理想气体模型，模拟Ar原子在特定温度下于二维平面上的热运动。

Ar原子的范德瓦尔斯半径为0.2 nm，因此可以粗略认为两个Ar原子在距离0.4 nm时发生了弹性碰撞。同时，当Ar原子距离容器壁距离小于0.2 nm时，我们认为其与容器壁发生了弹性碰撞。

这里设置如下初始条件：

所有分子的初始位置为随机值
所有分子的初始速率相同，速度方向为随机值

%matplotlib widgetimport numpy as np
from scipy.spatial.distance import pdist, squareform
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import matplotlib.path as path
from matplotlib.animation import FuncAnimationX, Y = 0, 1
class MDSimulation:def __init__(self, pos, vel, r, m):"""Initialize the simulation with identical, circular particles of radiusr and mass m. The n x 2 state arrays pos and vel hold the n particles'positions in their rows as (x_i, y_i) and (vx_i, vy_i)."""self.pos = np.asarray(pos, dtype=float)self.vel = np.asarray(vel, dtype=float)self.n = self.pos.shape[0]self.r = rself.m = mself.nsteps = 0def advance(self, dt):"""Advance the simulation by dt seconds."""self.nsteps += 1# Update the particles' positions according to their velocities.self.pos += self.vel * dt# Find indices for all unique collisions.dist = squareform(pdist(self.pos))iarr, jarr = np.where(dist < 2 * self.r)k = iarr < jarriarr, jarr = iarr[k], jarr[k]# For each collision, update the velocities of the particles involved.for i, j in zip(iarr, jarr):pos_i, vel_i = self.pos[i], self.vel[i]pos_j, vel_j = self.pos[j], self.vel[j]rel_pos, rel_vel = pos_i - pos_j, vel_i - vel_jr_rel = rel_pos @ rel_posv_rel = rel_vel @ rel_posv_rel = 2 * rel_pos * v_rel / r_rel - rel_velv_cm = (vel_i + vel_j) / 2self.vel[i] = v_cm - v_rel/2self.vel[j] = v_cm + v_rel/2# Bounce the particles off the walls where necessary, by reflecting# their velocity vectors.hit_left_wall = self.pos[:, X] < self.rhit_right_wall = self.pos[:, X] > 1 - self.rhit_bottom_wall = self.pos[:, Y] < self.rhit_top_wall = self.pos[:, Y] > 1 - self.rself.vel[hit_left_wall | hit_right_wall, X] *= -1self.vel[hit_bottom_wall | hit_top_wall, Y] *= -1# Number of particles.
n = 1000
# Scaling factor for distance, m-1. The box dimension is therefore 1/rscale.
rscale = 5.e6
# Use the van der Waals radius of Ar, about 0.2 nm.
r = 2e-10 * rscale
# Scale time by this factor, in s-1.
tscale = 1e9    # i.e. time will be measured in nanoseconds.
# Take the mean speed to be the root-mean-square velocity of Ar at 300 K.
sbar = 432 * rscale / tscale
# Time step in scaled time units.
FPS = 30
dt = 1/FPS
# Particle masses, scaled by some factor we're not using yet.
m = 1# Initialize the particles' positions randomly.
pos = np.random.random((n, 2))
# Initialize the particles velocities with random orientations and same magnitudes.
theta = np.random.random(n) * 2 * np.pi
s0 = sbar * np.ones(n)
vel = (s0 * np.array((np.cos(theta), np.sin(theta)))).Tsim = MDSimulation(pos, vel, r, m)# Set up the Figure and make some adjustments to improve its appearance.
DPI = 100
width, height = 800, 400
fig = plt.figure(figsize=(width/DPI, height/DPI), dpi=DPI)
fig.subplots_adjust(left=0, right=0.97)
sim_ax = fig.add_subplot(121, aspect='equal', autoscale_on=False)
sim_ax.set_xticks([])
sim_ax.set_yticks([])
# Make the box walls a bit more substantial.
for spine in sim_ax.spines.values():spine.set_linewidth(2)speed_ax = fig.add_subplot(122)
speed_ax.set_xlabel('Speed $v\,/m\,s^{-1}$')
speed_ax.set_ylabel('$f(v)$')particles, = sim_ax.plot([], [], 'ko')
particles_color, = sim_ax.plot([], [], 'ro')class Histogram:"""A class to draw a Matplotlib histogram as a collection of Patches."""def __init__(self, data, xmax, nbars, density=False):"""Initialize the histogram from the data and requested bins."""self.nbars = nbarsself.density = densityself.bins = np.linspace(0, xmax, nbars)self.hist, bins = np.histogram(data, self.bins, density=density)# Drawing the histogram with Matplotlib patches owes a lot to# https://matplotlib.org/3.1.1/gallery/animation/animated_histogram.html# Get the corners of the rectangles for the histogram.self.left = np.array(bins[:-1])self.right = np.array(bins[1:])self.bottom = np.zeros(len(self.left))self.top = self.bottom + self.histnrects = len(self.left)self.nverts = nrects * 5self.verts = np.zeros((self.nverts, 2))self.verts[0::5, 0] = self.leftself.verts[0::5, 1] = self.bottomself.verts[1::5, 0] = self.leftself.verts[1::5, 1] = self.topself.verts[2::5, 0] = self.rightself.verts[2::5, 1] = self.topself.verts[3::5, 0] = self.rightself.verts[3::5, 1] = self.bottomdef draw(self, ax):"""Draw the histogram by adding appropriate patches to Axes ax."""codes = np.ones(self.nverts, int) * path.Path.LINETOcodes[0::5] = path.Path.MOVETOcodes[4::5] = path.Path.CLOSEPOLYbarpath = path.Path(self.verts, codes)self.patch = patches.PathPatch(barpath, fc='tab:green', ec='k',lw=0.5, alpha=0.5)ax.add_patch(self.patch)def update(self, data):"""Update the rectangle vertices using a new histogram from data."""self.hist, bins = np.histogram(data, self.bins, density=self.density)self.top = self.bottom + self.histself.verts[1::5, 1] = self.topself.verts[2::5, 1] = self.topdef get_speeds(vel):"""Return the magnitude of the (n,2) array of velocities, vel."""return np.hypot(vel[:, X], vel[:, Y])def get_KE(speeds):"""Return the total kinetic energy of all particles in scaled units."""return 0.5 * sim.m * sum(speeds**2)speeds = get_speeds(sim.vel)
speed_hist = Histogram(speeds, 2.5 * sbar, 50, density=True)
speed_hist.draw(speed_ax)
speed_ax.set_xlim(speed_hist.left[0], speed_hist.right[-1])
# TODO don't hardcode the upper limit for the histogram speed axis.
ticks = np.linspace(0, 1000, 11, dtype=int)
speed_ax.set_xticks(ticks * rscale/tscale)
speed_ax.set_xticklabels([str(tick) for tick in ticks])
speed_ax.set_yticks([])fig.tight_layout()# The 2D Maxwell-Boltzmann equilibrium distribution of speeds.
mean_KE = get_KE(speeds) / n
a = sim.m / 2 / mean_KE
# Use a high-resolution grid of speed points so that the exact distribution
# looks smooth.
sgrid_hi = np.linspace(0, speed_hist.bins[-1], 200)
f = 2 * a * sgrid_hi * np.exp(-a * sgrid_hi**2)
mb_line, = speed_ax.plot(sgrid_hi, f, c='0.7')
# Maximum value of the 2D Maxwell-Boltzmann speed distribution.
fmax = np.sqrt(sim.m / mean_KE / np.e)
speed_ax.set_ylim(0, fmax)# For the distribution derived by averaging, take the abcissa speed points from
# the centre of the histogram bars.
sgrid = (speed_hist.bins[1:] + speed_hist.bins[:-1]) / 2
mb_est_line, = speed_ax.plot([], [], c='r')
mb_est = np.zeros(len(sgrid))# A text label indicating the time and step number for each animation frame.
xlabel, ylabel = sgrid[-1] / 2 + 0.1, 0.9 * fmax
label = speed_ax.text(xlabel, ylabel, '$t$ = {:.1f}s, step = {:d}'.format(0, 0))def init_anim():"""Initialize the animation"""particles.set_data([], [])particles_color.set_data([], [])return particles, particles_color, speed_hist.patch, mb_est_line, labeldef animate(i):"""Advance the animation by one step and update the frame."""global sim, verts, mb_est_line, mb_estsim.advance(dt)particles.set_data(sim.pos[:, X], sim.pos[:, Y])particles.set_markersize(1.0)particles_color.set_data(sim.pos[-5:, X], sim.pos[-5:, Y])particles_color.set_markersize(5.0)speeds = get_speeds(sim.vel)speed_hist.update(speeds)# Once the simulation has approached equilibrium a bit, start averaging# the speed distribution to indicate the approximation to the Maxwell-# Boltzmann distribution.if i >= IAV_START:mb_est += (speed_hist.hist - mb_est) / (i - IAV_START + 1)mb_est_line.set_data(sgrid, mb_est)label.set_text('$t$ = {:.1f} ns, step = {:d}'.format(i*dt, i))return particles, particles_color, speed_hist.patch, mb_est_line, label# Only start averaging the speed distribution after frame number IAV_ST.
IAV_START = 250
# Number of frames; set to None to run until explicitly quit.
frames = 301
anim = FuncAnimation(fig, animate, frames=frames, interval=10, blit=False,init_func=init_anim)anim.save('Maxwell-BoltzmannDistribution.gif', writer='imagemagick', fps=10)

效果如下

from IPython.display import Image
Image(filename ='Maxwell-BoltzmannDistribution.gif', width=800)

经过一定时间演化，速率分布变化逐步减小，达到平衡状态。其中灰线为二维情况下Maxwell-Boltzmann速率分布函数：
f(v)=mvkBTexp⁡(−mv22kBT)f(v) = \frac{mv}{k_\mathrm{B}T}\exp\left(-\frac{mv^2}{2k_\mathrm{B}T}\right)f(v)=kBTmvexp(−2kBTmv2)

1.4 Lennard-Jones potential

The potential energy of the 12-6 Lennard-Jones potential is given as

VLJ(r)=4ε[(σr)12−(σr)6]V_{\text{LJ}}(r)=4\varepsilon \left[\left({\frac {\sigma }{r}}\right)^{12}-\left({\frac {\sigma }{r}}\right)^{6}\right] VLJ(r)=4ε[(rσ)12−(rσ)6]

σ\sigmaσ is the radius where the potential is zero and is defined as the van der waals radius. ε\varepsilonε is the energy minimum of the interaction (see the figure below).

We don’t want to compute long range interactions as they will be negligible. We therefore apply a cutoff past rσ=2.5\frac{r}{\sigma}=2.5σr=2.5. However, as the particles move past the cutoff distance there will be a small jump in the energy, which is not realistic, so we shift the potential so that the potential goes to zero at rσ=2.5\frac{r}{\sigma}=2.5σr=2.5.

Below is a plot of the Lennard-Jones potential with the shifted potential shown for comparison.

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as npr = np.linspace(0.01,3.0, 500) # Make a radius vector
epsilon = 1 # Energy minimum
sigma = 1 # Distance to zero crossing point
E_LJ = 4*epsilon*((sigma/r)**12-(sigma/r)**6) # Lennard-Jones potentialplt.figure(figsize=[6,6])
plt.plot(r, E_LJ, 'r-', linewidth=1, label=r" $LJ\; pot$") # Red line is unshifted LJ# The cutoff and shifting value
Rcutoff = 2.5
phicutoff = 4.0/(Rcutoff**12)-4.0/(Rcutoff**6) # Shifts the potential so at the cutoff the potential goes to zeroE_LJ_shift = E_LJ - phicutoff # Subtract the value of the potential at r=2.5plt.plot(r[r < 2.5], E_LJ_shift[r < 2.5], 'b-', linewidth=1, label=r"$LJ\; pot\; shifted$") # Blue line is shifted#Plot formatting
plt.axhline(0, color='grey', linestyle='--', linewidth=2)
plt.axvline(1, color='grey', linestyle='--', linewidth=2)
plt.xlim([0.0, 3.0])
plt.ylim([-1, 1])
plt.ylabel(r"$E_{LJ}/\epsilon$", fontsize=20)
plt.xlabel(r"$r/\sigma$", fontsize=20)
plt.legend(frameon = False, fontsize=20)
plt.title(r"$Lennard-Jones \; potential$", fontsize=20)

1.5 Periodic boundary conditions

Periodic boundary conditions allow for an approximation of an infinitely sized system by simulating a simple unit cell. This is illustrated below. The black box is the only cell we simulate; the tiled images around it are there for illustration. The green particle moves past the top boundary of the unit cell and are moved to the bottom of the box with the same velocity (illustrated by the red dashed line). This boundary condition keeps the volume and number of particles constant in the simulation.

1.6 Reduced units

By choosing our units we can remove any constants and get a general behaviour for all gases. Mass, sigma, epsilon and the Boltzmann constant are set to equal one. Reduced coordinates are used for the other variables which are derived from the parameters set to one.
x∗=xσx^{*} = \frac{x}{\sigma} x∗=σx
v∗=vt∗σv^{*} = v\frac{t^{*}}{\sigma} v∗=vσt∗
t∗=t(ϵmσ2)1/2t^{*} = t\left(\frac{\epsilon}{m \sigma^{2}} \right)^{1/2} t∗=t(mσ2ϵ)1/2
E∗=EϵE^{*} = \frac{E}{\epsilon} E∗=ϵE
F∗=fσϵF^{*} = f\frac{\sigma}{\epsilon} F∗=fϵσ
P∗=Pσ3ϵP^{*} = P \frac{\sigma^{3}}{\epsilon}P∗=Pϵσ3
ρ∗=ρσdimensions\rho^{*} = \rho \sigma^{dimensions} ρ∗=ρσdimensions
T∗=TkbϵT^{*} = T \frac{k_{b}}{\epsilon} T∗=Tϵkb

Where kbk_{b}kb is Boltzmann’s constant.

This may seem complicated but it allows all of the equations to be written very simply in the program. This also gives us physical insight - for example, the reduced temperature is the ratio of the thermal energy kBTk_{B} TkBT to the energy of the intermolecular interactions ϵ\epsilonϵ.
Periodic boundary conditions

1.7 Calculating the forces

As mentioned above, the forces between particles can be calculated from the derivative/gradient of their potential energy. F=−1r∇E(r)\textbf{F}=-\frac{1}{r}\nabla E(\textbf{r})F=−r1∇E(r) (in spherical coordinates)

F=−1r∇ELJ(r)=−1rdELJdr=−24[2(σr)14−(σr)8]\textbf{F} = -\frac{1}{r}\nabla E_{LJ}(\textbf{r}) = -\frac{1}{r}\frac{d E_{LJ}}{d \textbf{r}} = -24\left[2\left(\frac{\sigma}{\textbf{r}}\right)^{14} - \left(\frac{\sigma}{\textbf{r}}\right)^{8}\right] F=−r1∇ELJ(r)=−r1drdELJ=−24[2(rσ)14−(rσ)8]

Periodic boundary conditions have to be considered when we compute the forces between particles because a particle near the boundary of the unit cell has to be able to feel the force from a particle on the other side of the unit cell. For example, the pink particle above will feel the force from the green particle, even though they are far from each other because they are near opposite boundaries.

In order to easily implement periodic boundary conditions, scaled box units are used so that the particle positions are always between -0.5 and 0.5. If the distance between the particles is greater than half the scaled box units, the interaction with the particle in the next box are consider

def Compute_Forces(pos, acc, ene_pot, epsilon, BoxSize, DIM, N):# Compute forces on positions using the Lennard-Jones potential# Uses double nested loop which is slow O(N^2) time unsuitable for large systemsSij = np.zeros(DIM) # Box scaled unitsRij = np.zeros(DIM) # Real space units#Set all variables to zeroene_pot = ene_pot*0.0acc = acc*0.0virial = 0.0# Loop over all pairs of particlesfor i in range(N-1):for j in range(i+1,N): #i+1 to N ensures we do not double countSij = pos[i,:]-pos[j,:] # Distance in box scaled unitsfor l in range(DIM): # Periodic interactionsif (np.abs(Sij[l])>0.5):Sij[l] = Sij[l] - np.copysign(1.0,Sij[l]) # If distance is greater than 0.5  (scaled units) then subtract 0.5 to find periodic interaction distance.Rij = BoxSize*Sij # Scale the box to the real units in this case reduced LJ unitsRsqij = np.dot(Rij,Rij) # Calculate the square of the distanceif(Rsqij < Rcutoff**2):# Calculate LJ potential inside cutoff# We calculate parts of the LJ potential at a time to improve the efficieny of the computation (most important for compiled code)rm2 = 1.0/Rsqij # 1/r^2rm6 = rm2**3.0 # 1/r^6rm12 = rm6**2.0 # 1/r^12phi = epsilon*(4.0*(rm12-rm6)-phicutoff) # 4[1/r^12 - 1/r^6] - phi(Rc) - we are using the shifted LJ potential# The following is dphi = -(1/r)(dV/dr)dphi = epsilon*24.0*rm2*(2.0*rm12-rm6) # 24[2/r^14 - 1/r^8]ene_pot[i] = ene_pot[i]+0.5*phi # Accumulate energyene_pot[j] = ene_pot[j]+0.5*phi # Accumulate energyvirial = virial + dphi*np.sqrt(Rsqij) # Virial is needed to calculate the pressureacc[i,:] = acc[i,:]+dphi*Sij # Accumulate forcesacc[j,:] = acc[j,:]-dphi*Sij # (Fji=-Fij)return acc, np.sum(ene_pot)/N, -virial/DIM # return the acceleration vector, potential energy and virial coefficient

1.8 Temperature

Temperature is a macroscopic quantity. Microscopically it is less well defineddue to the low number of particles. However, if we use the kinetic energy of the parameters we can calculate the temperature.
EK=12mv2E_{K}=\frac{1}{2} mv^{2} EK=21mv2kBT=23∑NEKk_{B} T = \frac{2}{3}\sum_{N}E_{K} kBT=32N∑EK

Where we sum over all $N $ atoms. We will use this in order to scale the velocities to maintain a constant temperature (remember we are using reduced units so kB=1k_{B}=1kB=1 and m=1m=1m=1).

The function below calculates the temperature given the velocity of the atoms/particles.

def Calculate_Temperature(vel,BoxSize,DIM,N):ene_kin = 0.0for i in range(N):real_vel = BoxSize*vel[i,:]ene_kin = ene_kin + 0.5*np.dot(real_vel,real_vel)ene_kin_aver = 1.0*ene_kin/Ntemperature = 2.0*ene_kin_aver/DIMreturn ene_kin_aver,temperature

1.9 Molecular dynamics program

The molecular dynamics program contains the instructions for the computer to use to move the particles/atoms through time. Most often this is written in Fortran and C; these compiled languages are orders of magnitude faster than Python, however a scripting language like Python is helpful to provide understanding about how molecular dynamics is implemented.

The main steps in a molecular dynamics simulation are:

Initialise the position of particles
Calculate the pairwise forces on the particles by calculating the gradient of the potential energy $ F = \nabla E(\textbf{r})=1/r\partial E(\textbf{r})/\partial r$
Compute the new positions by integrating the equation of motion (we will use the velocity Verlet algorithm)
Apply a thermostat to maintain the temperature at the set value (we will use the velocity scaling for temperature control)
Go back to step 2, recompute the forces and continue until the maximum number of steps

1.9.1 Initialising the particles

import numpy as npDIM = 2 # Dimensions
N = 32BoxSize = 8.0volume  = BoxSize**DIM
density = N / volume
print("volume = ", volume, " density = ", density)pos = np.random.rand(N, DIM)MassCentre = np.sum(pos,axis=0)/Nfor i in range(DIM):pos[:,i] = pos[:,i]-MassCentre[i]

1.9.2 Integrating the equations of motion

We will make use of the velocity Verlet integrator which integrates Newton’s equations of motion in 1D:
dxdt=vanddvdt=dF(x)m\frac{dx}{dt}=v\; and \; \frac{dv}{dt}=\frac{dF(x)}{m} dtdx=vanddtdv=mdF(x)

The velocity Verlet algorithm spilts the velocity update into two steps intially doing a half step then modifing the acceleration and then doing the second velocity update. Written in full, this gives:

Calculate x(t+Δt)=x(t)+v(t+12)Δtx(t+\Delta t) = x(t)+v\left(t+\frac{1}{2}\right)\Delta tx(t+Δt)=x(t)+v(t+21)Δt
Calculate $v\left(t+\frac{1}{2}\Delta t\right)= v(t)+\frac{1}{2}a(t)\Delta t $
Derive a(t+Δt)a(t+\Delta t)a(t+Δt) from the interaction potential using x(t+Δt)x(t+\Delta t)x(t+Δt)
Calculate v(t+Δt)=v(t+12)+12a(t+Δt)Δtv(t+\Delta t)=v\left(t+\frac{1}{2}\right)+\frac{1}{2}a(t+\Delta t) \Delta tv(t+Δt)=v(t+21)+21a(t+Δt)Δt

Between step 1 and 2 we rescale the velocities to maintain the temperature at the requested value.

The output is saved in the same folder as the ipython notebook as traj.xyz and can be opened by Avogadro or VMD.

# Main MD loop
def MD(pos, NSteps, deltat, TRequested, DumpFreq, epsilon, BoxSize, DIM):# Vectors to store parameter values at each stepN = np.size(pos[:,1])ene_kin_aver = np.ones(NSteps)ene_pot_aver = np.ones(NSteps)temperature = np.ones(NSteps)virial = np.ones(NSteps)pressure = np.ones(NSteps)ene_pot = np.ones(N)vel = (np.random.randn(N,DIM)-0.5)acc = (np.random.randn(N,DIM)-0.5)# Open file which we will save the outputs tof = open('traj.xyz', 'w')for k in range(NSteps):# Refold positions according to periodic boundary conditionsfor i in range(DIM):period = np.where(pos[:,i] > 0.5)pos[period,i]=pos[period,i]-1.0period = np.where(pos[:,i] < -0.5)pos[period,i]=pos[period,i]+1.0# r(t+dt) modify positions according to velocity and accelerationpos = pos + deltat*vel + 0.5*(deltat**2.0)*acc # Step 1# Calculate temperatureene_kin_aver[k],temperature[k] = Calculate_Temperature(vel,BoxSize,DIM,N)# Rescale velocities and take half stepchi = np.sqrt(TRequested/temperature[k])vel = chi*vel + 0.5*deltat*acc # v(t+dt/2) Step 2# Compute forces a(t+dt),ene_pot,virialacc, ene_pot_aver[k], virial[k] = Compute_Forces(pos,acc,ene_pot,epsilon,BoxSize,DIM,N) # Step 3# Complete the velocity step vel = vel + 0.5*deltat*acc # v(t+dt/2) Step 4# Calculate temperatureene_kin_aver[k],temperature[k] = Calculate_Temperature(vel,BoxSize,DIM,N)# Calculate pressurepressure[k]= density*temperature[k] + virial[k]/volume# Print output to file every DumpFreq number of stepsif(k%DumpFreq==0): # The % symbol is the modulus so if the Step is a whole multiple of DumpFreq then print the valuesf.write("%s\n" %(N)) # Write the number of particles to file# Write all of the quantities at this step to the filef.write("Energy %s, Temperature %.5f\n" %(ene_kin_aver[k]+ene_pot_aver[k],temperature[k]))for n in range(N): # Write the positions to filef.write("X"+" ")for l in range(DIM):f.write(str(pos[n][l]*BoxSize)+" ")f.write("\n")if(DIM==2):scatter.x = pos[:,0]*BoxSizescatter.y = pos[:,1]*BoxSizetime.sleep(0.1)#print(ene_kin_aver[k], ene_pot_aver[k], temperature[k], pressure[k]) f.close() # Close the filereturn ene_kin_aver, ene_pot_aver, temperature, pressure, pos

from bqplot import pyplot as plt
import time# Setting up the simulation
NSteps = 5000 # Number of steps
deltat = 0.0032 # Time step in reduced time units
TRequested = 0.5# #Reduced temperature
DumpFreq = 100 # Save the position to file every DumpFreq steps
epsilon = 1.0 # LJ parameter for the energy between particlesfigure = plt.figure(title='Molecular Dynamics')
figure.layout.width = '600px'
figure.layout.height = '600px'scatter = plt.scatter(pos[:,0]*BoxSize, pos[:,1]*BoxSize)plt.xlim(-0.5*BoxSize, 0.5*BoxSize)
plt.ylim(-0.5*BoxSize, 0.5*BoxSize)
plt.show()

ene_kin_aver, ene_pot_aver, temperature, pressure, pos = MD(pos,NSteps,deltat,TRequested,DumpFreq,epsilon,BoxSize,DIM)

import matplotlib.pyplot as plt# Plot all of the quantities
def plot():plt.figure(figsize=[7,12])plt.rc('xtick', labelsize=15) plt.rc('ytick', labelsize=15)plt.subplot(4, 1, 1)plt.plot(ene_kin_aver[1000:],'k-')plt.ylabel(r"$E_{K}$", fontsize=20)plt.subplot(4, 1, 2)plt.plot(ene_pot_aver[1000:],'k-')plt.ylabel(r"$E_{P}$", fontsize=20)plt.subplot(4, 1, 3)plt.plot(temperature[1000:],'k-')plt.ylabel(r"$T$", fontsize=20)plt.subplot(4, 1, 4)plt.plot(pressure[1000:],'k-')plt.ylabel(r"$P$", fontsize=20)plt.show()plot()

参考文献来自桑鸿乾老师的课件：科学计算和人工智能
特此鸣谢！！！

【Python分子动力学】相关推荐

使用Python代码实现一个简单的分子动力学模拟程序
1. 前言理解分子动力学模拟最好的方法是编写一个分子动力学程序,换句话说,教会计算机去理解并做执行,自己才算理解会了.因此本文将从常用于描述分子间的非键相互作用中的Lennard-Jones pot ...
AMBER:对单个复合物进行分子动力学模拟的python包(resp计算电荷及gpu加速版本)
保证cpu至少要有24个核或者自己改一下代码 #!/usr/bin/env python # -*- coding: utf-8 -*-# ============================= ...
python web游戏实例_python实现的简单文本类游戏实例
Python应用与实践 Python应用与实践目录 1. Python是什么? 1.1. Python语言 1.2. Python哲学 2. Python在工 ...
分子动力学模拟软件_基于GPU的分子动力学软件ACEMD的简介与安装
Acellera软件包括HTMD.ACEMD.AceCloud.Parameterize.AceFlow和ACEMD3模块. ACEMD简介 ACEMD是一款功能强大的生物分子动力学模拟软件包,该软件 ...
gromacs manual_GROMACS蛋白配体分子动力学模拟结果分析简要笔记
0. 引言本文以前文(https://zhuanlan.zhihu.com/p/149862369)为基础,对蛋白配体复合物分子模拟体系的结果进行一系列的粗浅分析,本文记述了简要的分析方法. 1 M ...
分子动力学模拟之周期性边界处理
技术背景周期性边界是分子动力学模拟中常用的一种技术手段,不仅可以完整的概述完整的分子体系的特性,在一部分场景下还可以提升计算的效率,从作用上来看更像是一类的近似模型(假设有一个原子逃出这个周期性边界 ...
值得推荐的分子动力学模拟入门书籍
来源:"分子动力学"公众号链接:https://mp.weixin.qq.com/s/a7SKL2Vbqv9MHaFueAg4tQ 入门阶段,首先你要知道你想做什么,最好是找个看 ...
分子动力学模拟之SETTLE约束算法
Python微信订餐小程序课程视频 https://edu.csdn.net/course/detail/36074 Python实战量化交易理财系统 https://edu.csdn.net/cou ...
分子动力学模拟gro格式转换为 car
分子模拟中,Materials Studio 跑分子动力学很慢,效率是gromacs的1/20.为了提高效率,可先在gromacs上运行分子动力学,获得平衡结构.再将gro文件转换为car文件并导入M ...
分子动力学开源分析软件MDAnalysis安装介绍及使用
MD Analysis安装和介绍一.MD Analysis 安装普通安装从源码安装: 二.MD Analysis使用读取文件 MD Analysis (以下简称mda) 是一款分子动力学后处理 ...

【Python分子动力学】