在周期分子动力学模拟中实现速度Verlet算法时的失控速度

Question

我正在尝试为具有周期性或反射边界条件的 2d 框中的 lennard-Jones 流体编写分子动力学模拟。模拟似乎在反射边界下运行良好，但由于某种原因，周期盒中的粒子在几百个积分步骤后开始获得额外的速度，最终没有粒子留在盒子里。

我知道失控速度可能是由于将时间步长设置得太小造成的。然而，当我使用相同的速度-verlet算法运行模拟时，我看不出这怎么可能是速度失控的原因，特别是当反射边界工作得很好时。我认为在这两种情况下，用于评估粒子力或粒子分离的方法可能存在错误，但我自己无法找到它。

我用来评估力的代码如下所示：

public static Vector3d LJForce(Particle3d a, Particle3d b,
        Vector3d particleSep) {
    //Define R (for simplicity)
    Vector3d R = particleSep; //(b-a)
    return new Vector3d(Vector3d.multVector3d(
            -24*EPSILON*(2*power12(SIGMA/R.getMag())/R.getMag()
                    -power6(SIGMA/R.getMag())/R.getMag()), R));
}

这应该赋予伦纳德-琼斯作用在每个粒子上的力。颗粒分离取决于边界条件。对于反射，它是：

@Override
public Vector3d getParticleSeparation(Particle3d a, Particle3d b) {
    return Particle3d.sepParticle3d(a, b);
}

与

@Override
public void checkBoundaries(Particle3d a)   {
    Vector3d position = a.getPosition();
    Vector3d velocity = a.getVelocity();
    /*
     * We want to check the x,y and z coordinates. If
     * they break the boundaries we need to reflect 
     * the position in that coordinate and invert the
     * velocity in that axis.
     */
    //X coordinate
    if(position.getX()<0.0) {
        position.setX(Math.abs(position.getX()));
        velocity.setX(-velocity.getX());    
    }
    else if(position.getX()>getWidth()) {
        double howFar = position.getX()-getWidth();
        position.setX(getWidth()-howFar);
        velocity.setX(-velocity.getX());
    }   
    else    {}
    //Y coordinate
    if(position.getY()<0.0) {
        position.setY(Math.abs(position.getY()));
        velocity.setY(-velocity.getY());    
    }
    else if(position.getY()>getWidth()) {
        double howFar = position.getY()-getWidth();
        position.setY(getWidth()-howFar);
        velocity.setY(-velocity.getY());
    }   
    else    {}
    //Z coordinate
    if(position.getZ()<0.0) {
        position.setZ(Math.abs(position.getZ()));
        velocity.setZ(-velocity.getZ());    
    }
    else if(position.getZ()>getWidth()) {
        double howFar = position.getZ()-getWidth();
        position.setZ(getWidth()-howFar);
        velocity.setZ(-velocity.getZ());
    }   
    else    {}
    a.setPosition(position);
    a.setVelocity(velocity);
}//End checkBoundaries(i)

将颗粒保存在盒子中。另一方面，对于周期性...

@Override
public Vector3d getParticleSeparation(Particle3d a,Particle3d b) {
    Vector3d sep = Particle3d.sepParticle3d(a, b);
    if(sep.getX()>=(double)getWidth()/2)    {
        sep.setX(sep.getX()-getWidth());
    }   else{}
    if(sep.getY()>=(double)getWidth()/2)    {
        sep.setY(sep.getY()-getWidth());
    }   else{}
    if(sep.getZ()>=(double)getWidth()/2)    {
        sep.setZ(sep.getZ()-getWidth());
    }   else{}
    return sep;
}

提供颗粒分离和

@Override
public void checkBoundaries(Particle3d a) {
    Vector3d position = new Vector3d();
    position.setX((getWidth()+a.getPosition().getX())%getWidth());
    position.setY((getWidth()+a.getPosition().getY())%getWidth());
    position.setZ((getWidth()+a.getPosition().getZ())%getWidth());
    a.setPosition(position);
}

检查边界。两种分离方法都调用

//Relative seperation
public static Vector3d sepParticle3d( Particle3d a, Particle3d b) {
    return new Vector3d(Vector3d.subVector3d(b.getPosition(),a.getPosition()));
}

这只是笛卡尔空间中两个粒子的向量分离。我的编码是否有错误或导致我的模拟因周期性边界条件而崩溃的其他原因？

谢谢

Answer 1

也许你已经想通了，但我认为错误在于周期边界的 checkBorders 代码。我可能是错的，但我认为你不能在浮点数/双精度上使用 % 运算符。 我认为代码应该是（另见 https://en.wikipedia.org/wiki/Periodic_boundary_conditions）

Vector3d op = a.getPosition(); //old position
double w = getWidth();
position.setX( op.getX() - w*((int)(2*op.getX()-1)) );
position.setX( op.getY() - w*((int)(2*op.getY()-1)) );
position.setX( op.getZ() - w*((int)(2*op.getZ()-1)) );