Simulating a More Interesting System
Now we consider a more interesting system consisting of many particles. We introduce a new member function entitled InitializeParticles? to add many particles to the simulation. We also change our previous local declaration of initial_conditions to a member variable (and rename it mIntialConditions), so that we can use it from all member functions.
#include <glotzilla++.h> class LennardJonesSimulation : public MdSimulation { public: LennardJonesSimulation() { InitializeParticles(); SetInteraction(new LennardJones); SetForceRoutine(new BruteForce); SetBoundaryConditions(new PeriodicBoundary(11)); SetMoveRoutine(new VerletMove); } private: void InitializeParticles() { ipstream particle_pipe("randbox -n=512 -r=-5,5,-5,5,-5,5 -m=1.0 -i=1000000"); vec_t x; while(particle_pipe >> x) { MdPointParticle *p = new MdPointParticle; p -> SetPosition(x); AddParticle(p); } particle_pipe.close(); } }; int main(int argc, char ** argv) { LennardJonesSimulation my_simulation; opstream my_pipe("vis3d -n=\"Lennard Jones System\""); for(int timesteps=0; timesteps<10000; timesteps++) { my_simulation +=5; my_simulation.PrintVisual(my_pipe); } my_pipe.close(); return 0; }
Now instead of entering particle positions manually, we read particle positions from a pipe.
ipstream particle_pipe("randbox -n=512 -r=-5,5,-5,5,-5,5 -m=1.0 -i=1000000");
Here, we open the Glotzilla "randbox" utility which creates random points in a box. The options specify that the program should try to create 512 points (-n=512) in a box that spans -5 to 5 in the x, y, and z directions (-r=-5,5,-5,5,-5,5). Note that this region is smaller than the box. This is so particles don't overlap on the boundary. The -m=1.0 option specifies that the points should be separated by at least 1.0. The i=1000000 option specifies that the program should give up trying to add points after 1000000 overlaps.
When randbox outputs points to stdout, they captured and put into the system using the following commands:
vec_t x;
while(particle_pipe >> x)
{
MdPointParticle *p = new MdPointParticle;
p -> SetPosition(x);
mInitialConditions -> AddParticle(p);
}