Runaway Velocities when implementing Velocity Verlet Algorithm in Periodic Molecular Dynamics Simulation

Question

I'm trying to write a molecular dynamics simulation for a lennard-jones fluid in a 2d box with either periodic or reflective boundary conditions. The simulation seems to run well with reflective boundaries but for some reason the particles in the periodic box start to gain extra velocity after a few hundred integration steps and eventually none of the particles remain in the box.

I understand that runaway velocities can be caused by setting the time-step too small. However as I run the simulation with the same velocity-verlet algorithm I can't see how this could be the reason for the runaway velocities, especially as reflective boundaries work so well. I think there might be a mistake in the methods used to evaluate the forces on the particles or the particle separations in the two cases but I can't locate it myself.

The code I'm using to evaluate the forces looks like this:

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));
}

This is supposed to give the lennard-jones force acting on each particle. The particle separation is determined depending on the boundary conditions. For reflective it is:

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

which is used in conjunction with

@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)

to keep the particle in the box. For periodic on the other hand...

@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;
}

gives the particle separation and

@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);
}

checks the boundaries. Both of the separation methods call

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

which is just the vector separation of the two particles in Cartesian space. Is there a mistake in my coding or some other reason why my simulation breaks down with periodic boundary conditions?

Thanks

Answer 1

Maybe you figured it out already, but I think the mistake is in the checkBoundaries code for the periodic boundaries. I might be wrong, but I don't think you can use the % operator on floats/doubles. I think the code should be (see also 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)) );