If anybody stumbles upon this and has a similar problem/need, here's some solutions I've gathered. All pure Java.

• Traer Physics: created for use in Processing, is, in fact, pure Java. Easy to understand and extend. Implements Euler, Modified Euler and Runge-Kutta 4 integrators. I've written a Verlet Velocity integrator class, that's almost as fast as Euler but more precise (will put on github). A good middle way between Euler and RK4.
• Toxilibs: a more advanced collection of utilities for physics, collision detection etc. Will probably end up using this if I stop writing my own. It uses Verlet Integration.
• And of course, last and (for now) least, my own dirty little physics engine. I'm trying to write the simplest and fastest version without caring too much about readability/code reuse. If you are trying to learn these things like, it might be useful to study my code. Or not.

You can do this yourself pretty easily. I imagine your simulation of forces can be very simple and straightforward. To me it sounds like the collision detection would be the hardest part, but luckily for just points you can do this very easily.

You can use for a broadphase a hashed grid. The idea is extremely simple to implement and very fast. You can look up source code on how to construct a hashed grid in Ericson's Orange Book on Real Time Collision Detection.

The idea is to take the floating point components (x, y, z) of a point and then hash them into an integer index. The grid isn't stored in memory. Here's an example of finding out which grid cell an arbitrary point is in:

#define GRID_CELL_SIZE = 5.0
Vec3 particlePosition = particle.position;

int x = particlePosition.x / GRID_CELL_SIZE;
int y = particlePosition.y / GRID_CELL_SIZE;
int z = particlePosition.z / GRID_CELL_SIZE;

int hash = CoordinateHash( x, y, z );

// Example of hashing a coordinate
int CoordinateHash( int x, int y, int z )
{
    const int someLargePrimeNumber0;
    const int someLargePrimeNumber1;
    const int someLargePrimeNumber2;

    return x * someLargePrimeNumber0 + y * someLargePrimeNumber1 + z * someLargePrimeNumber2;
}

This is used to place the point into a hash table bucket. The hash collisions can be resolved through chaining.

To check for collision in the broadphase, hash a point and check it's chain for near-proximity with any other points in the chain (or check collision for a pair of points in any way you desire). Please note it doesn't matter if two different grid cells map to the same bucket, since you'd be doing a proximity check on a collision chain anyways.

Yep, in fact I've implemented a small Verlet Velocity integration in a few hours (mostly spent reading around and refreshing my physics).

For collisions, I was actually thinking on using (or implementing) a KD-Tree allowing range searches. As far as I understand, for points that's about the best. I haven't yet looked into it though: is the approach you mention comparable or totally unrelated?