如何在吸引粒子的模拟中添加硬球斥力？

Question

我模拟了一个 2D 粒子系统，它们相互吸引。吸引力的强度取决于粒子之间的距离。边界条件和相互作用是周期性的。由于吸引力，粒子相互靠近并聚集成一个圆圈。

我想添加硬球斥力，这样，每当两个或多个粒子聚集在同一位置时，我希望它们在连接它们中心的线上分开，直到它们不重叠。我该怎么做？

当存在吸引相互作用时添加硬球的情况比通常情况下更难，因为在某些情况下，同一位置可能有 4 个或更多粒子。

这是我的代码：

#include <iostream>
#include <math.h>
#include <vector>
#include <array>
#include <list>
#include <random>
#include <functional>
#include <fstream>
#include <string>
#include <sstream>
#include <algorithm>
#include <chrono>
#include <set>

using namespace std;

    std::minstd_rand gen(std::random_device{}());
    std::uniform_real_distribution<double> unirnd(0, 1);

double PBCtwo(double pos, double L)
    {
    if(pos > L / 2.0)
        return pos-L;
    else if (pos < -L /2.0)
        return L + pos;
    else
        return pos;
    }

// main function
int main()
{   long c = 0;
    int N=4000; 
    double rho, v0, tr,xr,l0, eta,dt, x[N],y[N],L=pow(N / rho , 0.5),l0_two = l0 * l0;
                rho = 2;
                v0 =300;eta = 1;dt = 0.0001;l0 = 1; c_prod = 500;c_display = 100;tr = -0.4;        
    // write initial configuration to the file
    ofstream configFile;
    configFile.open ("Initial Configuration.txt");
    configFile << to_string(N) << "\n";
    configFile << to_string(L) << "\n";
    for (int i = 0; i < N; i++)
    {   x[i] = unirnd(gen) * L;
        y[i] = unirnd(gen) * L;
        configFile << to_string(x[i]) << "\t" << to_string(y[i]) <<   "\n";
    }
    configFile.close();

    while (c < c_prod)
        {
        double dx[N], dy[N];
        c++;
        for(int i = 0; i < N; i++)
           {
            dx[i] = 0;
            dy[i] = 0;
            double S_try = 0.0, S_trx = 0.0;
            for(int j = 0; j < N; j++)
              {
                if (j==i) continue;
                   double delta_x = x[i]-x[j],  
                          delta_y = y[i]-y[j];
                       double r_x_ij = PBCtwo(delta_x,L),
                              r_y_ij = PBCtwo(delta_y,L),
                              r_ij_square = r_x_ij * r_x_ij + r_y_ij * r_y_ij;
                       if (r_ij_square > l0_two)
                          { 
                          double r_ij = sqrt(r_ij_square);
                          r_x_ij/= r_ij; 
                          r_y_ij/= r_ij;
                          double S_tr = 1 /r_ij_square;
                          S_trx += r_x_ij * S_tr;
                          S_try += r_y_ij * S_tr;
                          }
                }
            dx[i] += tr * S_trx;  
            dy[i] +=  tr * S_try;
            }
        for(int i = 0; i < N; i++)
            {
            x[i]+=  dt * dx[i];
            y[i]+=  dt * dy[i];
            if (x[i] > L){
              x[i]-= L;} 
            else if( x[i] < 0) {
            x[i]+= L;}
            if (y[i] > L){
             y[i]-= L;} 
            else if( y[i] < 0){
                y[i]+= L;}
            }
        }
    ofstream finalConfigFile;
    finalConfigFile.open ("Final Configuration.txt");
    finalConfigFile << to_string(N) << "\n";
    finalConfigFile << to_string(L) << "\n";
    for (int i = 0; i < N; i++)
        {
        finalConfigFile << to_string(x[i]) << "\t" << to_string(y[i]) <<"\n";   
        }
    finalConfigFile.close();
    return 0;
}

Answer 1

我假设你想模拟一种具有潜力的粒子气体，它具有吸引和核心排斥的部分。

多粒子系统的分子动力学模拟是一个发达的科学领域。这种模拟的难点在于最初以不重叠的方式对粒子进行采样。 Metropolis et al 在经典论文中解决了这个问题。但是无论如何你都忽略了这部分。还有很多教科书和论文解释了如何用硬球进行分子动力学模拟，比如this或this或1959 paper, where they first used Molecular Dynamics。

如果您的项目是像在 Windows 屏幕保护程序中一样显示碰撞的漂亮球体，只需让您的核心潜力更柔和，减少您的时间步长，它就会像魅力一样工作。如果您的目标是科学严谨，请阅读那些论文并且不要忘记控制能量守恒。可以添加软排斥，例如，通过替换

double S_tr = 1 /r_ij_square;

通过

double S_tr = 1.0 / r_ij_square - 1.0 / (r_ij_square * r_ij_square);

代码格式、结构和变量命名还有很多不足之处。我建议尝试clang-format 工具。在结构上，周期性边界条件有很多重复。

借助排斥势，您可能会得到一个经历相变的系统，这是一个令人兴奋的研究对象。祝你好运！

更新：

请注意，以下代码存在错误：

double rho, v0, tr, xr, l0, eta, dt, x[N], y[N];
double L = pow(N / rho , 0.5), l0_two = l0 * l0;
// ...
l0 = 1;

你首先计算 l0_two = l0 * l0（永远不要以 l 开头变量名，它很容易与 1 混淆，sqrt(x) 应该优先于 pow(x, 0.5)）然后分配 l0 = 1。因此 l0_two 是未定义的，可能被分配为 0。这可能会导致您的部分问题。

修复该错误后，您可以选择两种解决方案。如果粒子足够接近，第一个将是上述软排斥势。第二个是完全按照paper with event-based MD cited above 做的。

Answer 2

模拟这个硬球最常见的方法是，考虑到每个粒子都有一个半径，如果两个粒子之间的距离小于半径之和（它们在空间中重叠），则将它们的中心与一个（虚拟）连接起来春天会把它们带走。

斥力的公式类似于 F=-k*(l-l_0) 其中 k 是修改接触硬度的经验常数，l_0 是两个粒子半径之和，l 是实际两个中心之间的距离。仅当两个粒子重叠时才考虑此力！

为此，您必须时刻测试每对粒子是否重叠。如果是这样，请计算弹簧的力并将其考虑在内。