计算科学 - 讨论多原子的朗之万动力学模拟的能谱 - 吾爱随笔录

讨论多原子的朗之万动力学模拟的能谱

计算科学计算物理学模拟分子动力学

2021-12-02 09:12:31

更新

我已经在 3D 中编写了多粒子 MD 模拟。它基于朗之万动力学，具有随机脉冲和耗散。我认为该程序现在可以正常工作了吗？我附上了动能图、势能图和总能量图。我相信在没有朗之万动力学的情况下，总能量应该在耗散关闭时守恒。但是当耗散打开时，人们不应该再指望总能量守恒了。在这一点上，我想知道我的结果是否合理。我已经包含了一些更重要的代码片段（如果需要，请询问更多）

物理参数：

//V0 is the potential depth, r0 is the effective radius for the potential
double const m = 1., V0 = 1., r0 = 1., boxLength = 15., kT = 0.5;
int const n_atoms = 12;
double gamma = 0.5;   //damping coeff.
double var = 2.*gamma*m*kT*dt;  //variance of the gaussian distributed I
double c = (2.-gamma*dt)/(2.+gamma*dt);    //due to damping
double dt = 0.005;   //time step
int N = (int)(1000 / dt);   //simulation time

初始化系统

void init(struct Atom system[n_atoms]){
    double rc = 0.9;
    double xx = 0., yy = boxLength, zz = boxLength;
    int atoms_in_row = (int)(boxLength / rc);
    int n = 0;

    for(int k = 0; k < (int)boxLength; k++) {
        for(int i = 0; i <= atoms_in_row; i++) {
            for(int j = 0; j <= atoms_in_row; j++) {
                if(n >= n_atoms) break;
                system[n].x = xx;
                system[n].y = yy;
                system[n].z = zz;
                system[n].vx = system[n].vy = system[n].vz = 0.;
                system[n].ax = system[n].ay = system[n].az = 0.;
                xx += rc;
                n++;
            }
            yy -= rc;
            xx = 0.;
        }
        yy = 0.;
        zz -= rc;
    }
}

这是分别关闭和打开耗散/随机脉冲的能量图。第一个仅显示了非 LD 情况下两个原子的能量图。接下来的两个用于 12 个原子。特别是，我认为非 LD 情况的结果很好，因为总能量或多或少是守恒的。也许结果对 LD 案例也有好处？我不太确定这一点。

2个回答

根据您的评论，我认为您应该花更多时间调试代码。您说过有时程序会陷入无限循环：这是您需要解决的第一件事。

我建议您阅读三本关于分子动力学方法的非常清晰且写得很好的书：

1 D. Frenkel & B. Smit，了解分子模拟：从算法到应用，第 2 版。

[2] Allen & Tildesley，液体的计算机模拟。

[3] DC Rapaport，分子动力学模拟的艺术

它们带来了实现的例程，您可以在其中与您的代码进行比较，并验证您是否将方法很好地转换为编程。按照第一个作为 Molecular Dynamics 的主要参考，如果您可能有一些 Frenkel 中未阐明的疑问，请查看其他的。

Lennard-Jones 分子动力学参数化[1]：

$U(r) = 4 \epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right]$

tstep = 0.01; // time step
temperature = 0.728; 
epsilon = 1.0; // potential depth
sigma = 1.0; // effective radius
cut_off = 2.5*sigma; // LJ cut-off
min_distance = 0.87; // minimum initial distance between atoms

基于粗略查看您的代码，我认为您可能需要更改以下行calculateForce：

        system[i].ax += f * dx / d;
        system[i].ay += f * dy / d;
        system[i].az += f * dz / d;
        system[j].ax += -system[i].ax;
        system[j].ay += -system[i].ay;
        system[j].az += -system[i].az;

并用以下效果替换它们：

        system[i].ax += f * dx / d;
        system[i].ay += f * dy / d;
        system[i].az += f * dz / d;
        system[j].ax -= f * dx / d;
        system[j].ay -= f * dy / d;
        system[j].az -= f * dz / d;

不管你当前的问题是什么，一个好的调试策略是关闭朗之万动力学引入的波动和耗散（即，简单地使用 Verlet），看看你是否节约能源。如果答案是否定的，则意味着您需要修复一个错误。如果您不这样做，随机动态可能会使您更难发现代码中的问题。

其它你可能感兴趣的问题

上一篇模拟带障碍物的斯托克斯流下一篇是否可以用有限差分 WENO 方法求解多分量欧拉方程？