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

Question

我模拟了一个相互吸引的二维粒子系统。吸引力的强度取决于粒子之间的距离。边界条件和相互作用是周期性的。因为有吸引力，粒子会相互靠近，聚集成一个圆圈。

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

当存在吸引相互作用时添加硬球的情况比通常情况更难，因为在某些情况下可能有 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

最常见的模拟硬球体的方法是，考虑到每个粒子都有一个半径，如果两个粒子之间的距离小于半径之和（它们在space中重叠）link 他们的中心有一个（虚拟的）spring 会把他们带走。

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

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

Answer 2

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

多粒子系统的分子动力学模拟是一个发达的科学领域。这种模拟的困难部分是最初以不重叠的方式对粒子进行采样。这在经典论文中由Metropolis et al. But you neglect this part anyway. There are also many textbooks and papers explaining, how to perform molecular dynamics simulations with hard spheres, like this or this or the 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.

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

How to add hard-sphere repulsion in my simulation of attractive particles?

c++

simulation

c++11