Homework 0: Molecular Dynamics 分子动力学模拟

import numpy as np
import taichi as ti

ti.init(arch=ti.gpu)

# dimension of simulation system
DIM = 3
WINDOW_SIZE = 1024
rgb2hex = lambda x: x[0].astype(np.int) * 65536 + x[1].astype(np.int) * 256 + x[2].astype(np.int)

@ti.data_oriented
class MolecularDynamics:
    # cutoff radius, interaction is not calculated
    # if two particles are farther than cutoff
    rcut = 2.5
    ecut = 4. * (1. / rcut ** 12 - 1. / rcut ** 6)
    # thermostat "damping" coefficients
    Q1 = 5
    Q2 = 5

    '''
    Initializes the object, set up python scope
    variables and taichi vectors.
    '''
    def __init__(self, density, temperature, boxlength, dt=5e-3):
        # python scope variables
        self.rho = density
        self.boxlength = boxlength
        self.dt = dt
        self.n_particles = int(density * boxlength ** DIM)
        # ti variables
        self.position = ti.Vector(DIM, dt=ti.f32, shape=self.n_particles)
        self.velocity = ti.Vector(DIM, dt=ti.f32, shape=self.n_particles)
        self.force = ti.Vector(DIM, dt=ti.f32, shape=self.n_particles) 
        # kinetic energy
        self.ek = ti.var(dt=ti.f32, shape=())
        # potential energy
        self.ep = ti.var(dt=ti.f32, shape=())
        # temperature
        self.temp = ti.var(dt=ti.f32, shape=())
        # initialize system
        self.init_thermostat()
        self.init_pos()
        self.init_velocity()
        self.temp.fill(temperature)
        # spawns GUI
        self.gui = ti.GUI("MD", res=WINDOW_SIZE)
    
    def set_temp(self, temperature):
        self.temp.fill(temperature)

    '''
    Initializes the simulation system. Place particles on a regular grid
    and randomize their velocities according to the temperature. Also 
    initializes the thermostat.
    '''
    def init_pos(self):
        n_pow = int(self.n_particles ** (1. / DIM))
        # n_axes = [nx, ny, ...] is the number of particles along each axis to be placed.
        n_axes = np.array([n_pow] * DIM)
        for i in range(DIM):
            if n_pow ** (DIM - i) * (n_pow + 1) * i < self.n_particles:
                n_axes[i] += 1
        disp = self.boxlength / n_axes
        coords_1d = [d * (0.5 + np.arange(n)) for d, n in zip(disp, n_axes)]
        self.position.from_numpy(
            np.stack(np.meshgrid(*coords_1d)).reshape(DIM, -1)[:, :self.n_particles].T)

    def init_velocity(self):
        vs = np.random.random((self.n_particles, DIM)) - 0.5
        vcm = np.mean(vs, axis=0).reshape((1, DIM))
        vs -= vcm
        vs *= np.sqrt(DIM * self.temp[None] * self.n_particles / np.sum(vs ** 2))
        self.velocity.from_numpy(vs)

    def init_thermostat(self):
        self.xi = ti.var(dt=ti.f32, shape=2)
        self.vxi = ti.var(dt=ti.f32, shape=2)
        self.xi.fill(0)
        self.vxi.fill(0)

    '''
    Calculates distance with periodic boundary conditions
    and wraps a particle into the simulation box.
    '''
    @ti.func
    def calc_distance(self, x1, x2):
        dist = ti.Vector([0.0] * DIM)
        for i in ti.static(range(DIM)):
            dist[i] = x1[i] - x2[i]
            if dist[i] <= -0.5 * self.boxlength:
                dist[i] += self.boxlength
            elif dist[i] > 0.5 * self.boxlength:
                dist[i] -= self.boxlength
        return dist

    @ti.func
    def wrap(self, x):
        for i in ti.static(range(DIM)):
            if x[i] <= 0:
                x[i] += self.boxlength
            elif x[i] > self.boxlength:
                x[i] -= self.boxlength
        return x

    '''
    Calculates the interaction energy between two particles.
    '''
    @ti.func
    def interaction_energy(self, r2):
        return 4. * (1. / r2 ** 6 - 1. / r2 ** 3) - self.ecut

    '''
    Calculates magnitude of the force between two particles 
    multipiled by their distance.
    '''
    @ti.func
    def force_times_dist(self, r2):
        return 24. * (2. / r2 ** 6 - 1. / r2 ** 3)

    '''
    Substeps to integrate the thermostat.
    '''
    @ti.func
    def step_vxi1(self, G1):
        self.vxi[0] *= ti.exp(-self.vxi[1] * self.dt * 0.125)
        self.vxi[0] += G1 * self.dt * 0.25
        self.vxi[0] *= ti.exp(-self.vxi[1] * self.dt * 0.125)

    @ti.func
    def step_vxi2(self, G2):
        self.vxi[1] = self.vxi[1] + G2 * self.dt * 0.25

    @ti.func
    def step_xi(self):
        self.xi[0] += self.vxi[0] * self.dt * 0.5
        self.xi[1] += self.vxi[1] * self.dt * 0.5
    
    '''
    Integrate the thermostat by half a time step.
    '''
    @ti.kernel
    def integrate_thermostat(self):
        G2 = self.Q1 * self.vxi[0] ** 2 - self.temp[None]
        self.step_vxi2(G2)
        G1 = (2 * self.ek[None] - 3 * self.n_particles * self.temp[None]) / self.Q1
        self.step_vxi1(G1)
        self.step_xi()
        s = ti.exp(-self.vxi[0] * self.dt * 0.5)
        self.ek[None] *= s * s
        G1 = (2 * self.ek[None] - 3 * self.n_particles * self.temp[None]) / self.Q1
        self.step_vxi1(G1)
        G2 = (self.Q1 * self.vxi[0] ** 2 - self.temp[None]) / self.Q2
        self.step_vxi2(G2)
        for i in self.velocity:
            self.velocity[i] = self.velocity[i] * s

    '''
    Integrate the motion of particles. Use Newton'w law of motion and 
    verlet integration scheme. Also calculates the kinetic and potential energies.
    '''
    @ti.kernel
    def integrate_motion(self):
        self.ek[None] = 0.0
        self.ep[None] = 0.0
        for i in self.position:
            self.position[i] = self.wrap(self.position[i] \
                        + self.velocity[i] * self.dt * 0.5)
            self.force[i].fill(0)
        for i, j in ti.ndrange(self.n_particles, self.n_particles):
            if i < j:
                d = self.calc_distance(self.position[i], self.position[j])
                r2 = (d ** 2).sum()
                if r2 < self.rcut ** 2:
                    rf_norm = self.force_times_dist(r2)
                    # += performs atomic add
                    self.force[i] += rf_norm * d / r2
                    self.force[j] -= rf_norm * d / r2
                    self.ep[None] += self.interaction_energy(r2)

        for i in self.position:
            self.velocity[i] = self.velocity[i] + self.force[i] * self.dt
            self.ek[None] += (self.velocity[i] ** 2).sum() / 2
            self.position[i] = self.wrap(self.position[i] \
                        + self.velocity[i] * self.dt * 0.5)

    '''
    Renders a simulation frame to the GUI.
    '''
    def render(self, savefile=None):
        camera = self.boxlength
        radius = 8
        bg = np.array([17, 47, 65])
        circ = np.array([122, 200, 225])
        self.gui.clear(rgb2hex(bg))
        while self.gui.get_event(ti.GUI.PRESS):
            if self.gui.event.key == ti.GUI.ESCAPE:
                exit()
            if self.gui.event.key == ti.GUI.UP:
                self.temp[None] = round(self.temp[None] + 0.1, 1)
            if self.gui.event.key == ti.GUI.DOWN:
                if self.temp[None] <= 0.15:
                    self.temp[None] /= 2
                else:
                    self.temp[None] = round(self.temp[None] - 0.1, 1)
        xy = self.position.to_numpy()        
        z = xy[:, 2]
        xy = xy[:, :2]
        z_order = np.argsort(z)[::-1]
        sizes = radius * camera / (camera + z)
        z = (1 - z / np.max(z))
        colors = rgb2hex(np.outer(bg, 1 - z) + np.outer(circ, z))
        t_actual = 2 * self.ek[None] / (self.n_particles * DIM)
        color_t = 0xf56060 if abs(t_actual - self.temp[None]) > 0.02 else 0x74e662
        self.gui.circles(xy / self.boxlength, radius=sizes[z_order], color=colors[z_order])
        self.gui.text("T_set = %.3g" % self.temp[None], (0.05, 0.2), font_size=36)
        self.gui.text("T_actual = %.3g" % t_actual, (0.05, 0.15), font_size=36, color=color_t)
        self.gui.text("Internal energy = %.3f" % (self.ek[None] + self.ep[None]), (0.05, 0.1), font_size=36)
        self.gui.show(savefile)
    
    '''
    Runs the simulation.
    '''
    def run(self, time, irender=10, save=False):
        nframe = int(time / self.dt)
        for i in range(nframe):
            self.integrate_thermostat()
            self.integrate_motion()
            self.integrate_thermostat()
            if irender > 0 and i > 0 and i % irender == 0:
                if save:
                    self.render("frame%i.png" % i) 
                else:
                    self.render()

if __name__ == "__main__":
    # number of particles
    n = 4000
    density = 0.1
    temperature = 0.5
    boxlength = (n / density) ** (1 / DIM)
    md = MolecularDynamics(density, temperature, boxlength)
    print(md.n_particles, boxlength)
    # initially run a few steps at a high temperature
    # to randomize the structure
    md.set_temp(1.5)
    md.run(md.dt * 5000, irender=-1)
    md.set_temp(temperature)
    # Press UP to increase the temperature and press DOWN to decrease
    md.run(1e8)

分子动力学模拟(MD)通过牛顿运动定律计算微观体系中原子和分子的运动轨迹。和弹簧质点体系一样,都是求解运动方程对时间的积分。分子动力学一般使用显式的Verlet积分法, 除了计算粒子的速度和位置以外,为了在恒定温度下进行模拟还需要引入一组附加量,其作用是改变粒子的速度(Nose–Hoover thermostat).

9 Likes

单原子分子(如稀有气体)之间的相互作用可以用Lennard–Jones potential描述,在不同温度下可以模拟固体、液体、气体三种状态:


固态:固体分子的排列相对规律,在平衡位置附近振动。只有极少数的游离气态分子在体系中共存。

4 Likes


液态:液体分子的运动范围更大,故液体具有流动性,体系存在气液共存,气体区域的密度明显小于液体。

3 Likes


气态:气体分子能够在自由运动,整个体系的密度是均匀的。

:joy: :joy: :joy:一层楼只能发一个图也是十分尴尬了hhhh

6 Likes

Hi, 我有一个问题想请教。如果直接用 Lennard–Jones 力做流体模拟会怎样?
SPH方法利用NS方程,需要计算压强p然后Δp,如果直接用 LJ力不就可以跳过这一步了吗?这样计算简单,而且它某种程度上还能模拟出表面张力。
我知道LJ力是分子层面的表述,NS方程是宏观层面的表述。如果我们把particle当成一个水滴,把NS方程中的Δp换成LJ力,会有什么问题?

在SPH的原始版本里,lennard-jones力被用来当做penalty force来处理边界条件,防止粒子穿越边界,这种方法会造成粒子间的压力差别大,分布不平滑,容易不稳定,计算时需要选择很小的积分时间步长,我的理解是直接用在计算压力里也会有相似的问题。

1 Like