Initialize System on a Lattice

This tutorial creates a custom class called LennardJonesSimulationLattice that inherits from the MdSimulation class. This code uses a pipe stream to access the "mkcrystal" routine to initialize particles on a simple cubic lattice. Velocities are also randomly assigned using the random number generation class. This tutorial builds on the topics in Minimal Simulation OOP, so please read this first.

The following code pertains to a single component Lennard-Jones, NVE molecular dynamics simulation.

• Verified10-4-07; code revision 390
• Created 5-12-07; code revision 185

Inheriting from MdSimulation

We will start with the code from Minimal Simulation OOP:

#include <glotzilla++.h>

class LennardJonesSimulationLattice : public MdSimulation
{
   public:
      LennardJonesSimulationLattice()
      {                                  
                AddParticle(new MdPointParticle(0,0,0));
                AddParticle(new MdPointParticle(1,1,1));
                SetInteraction(new LennardJones);
                SetForceRoutine(new BruteForce);
                SetBoundaryConditions(new PeriodicBoundary(10));
                SetIntegrationScheme(new VelocityVerletIntegrator);                             
      }
};

int main(int argc, char ** argv) 
{
   LennardJonesSimulationLattice my_simulation;

   for(int timesteps=0; timesteps<10000; timesteps++)
      my_simulation ++;

   return 0;
}

We will focus on the constructor of our inherited class and change how particles are initialized.

Initialize Particles Function

We will remove the AddParticle? call from the constructor and replace it with a new custom function called "InitializeParticles".

class LennardJonesSimulationLattice : public MdSimulation
{
   public:
      LennardJonesSimulationLattice()
      {                                  
                InitializeParticles();
                SetInteraction(new LennardJones);
                SetForceRoutine(new BruteForce);
                SetBoundaryConditions(new PeriodicBoundary(10));
                SetIntegrationScheme(new VelocityVerletIntegrator);                             
      }
      void InitializeParticles()
      {

      }
};

mkcrystal

Within the "InitializeParticles" routine, we could write our own code to loop over the box and place particles on a lattice. However, Glotzilla provides routines to accomplish within the DataGeneration? Library; specifically we will be using the "mkcrystal" application. The easiest way to access this functionality is by use of a pipe stream; we will essentially be running a binary application and capturing the output.

The application "mkcrystal" has the following usage (accessible by entering "mkcrystal -h" in your terminal):

usage:

-h      --help                  Prints this help
-n      --numpart=INTEGER       Number of particles
-d      --density=DECIMAL       Number density of the system
-s      --shape=SHAPE           Shape of the crystal (default is cube)
-S      --sc                    Simple cubic
-F      --fcc                   Face centered cubic
-H      --hcp                   Hexagonal close packed

Pipestream

To access the output from "mkcrystal" we will utilize a pipestream. The following excerpt will call mkcrystal and create a 200 points arranged in a simple cubic lattice with a number density of 0.85 (this dictates the spacing between the points).

void InitializeParticles()
{
        ipstream particle_pipe("mkcrystal -n=200 -d=0.85 -S");
 
}

We will then need to loop over the output from mkcrystal that is stored in "particle_pipe". To do this we will will use a while loop where at each stage we store the output from each line into the vector "x" (mkcrystal outputs <x,y,z> coordinates per point, per line).

void InitializeParticles()
{
        ipstream particle_pipe("mkcrystal -n=200 -d=0.85 -S");
        gvector x;         
        while(particle_pipe >> x)
        {
                MdPointParticle *p = new MdPointParticle;                  
                p -> SetPosition(x);
                AddParticle(p);
        }
        particle_pipe.close();
}

Let us exam this code. We first create an instance of MdPointParticle "p". We then call the SetPosition? function of MdPointParticle passing it the vector "x" which contains the output from "mkcrystal"; SetPosition? can take either a vector or 3 float values as input. Finally we need to add this particle to the system by way of the AddParticle? function. Again, AddParticle? is inherited from MdSimulation. Recall how we previously operated on AddParticle?:

        AddParticle(new MdPointParticle(0,0,0));
        AddParticle(new MdPointParticle(1,1,1));

Setting Particle Velocity

Within the InitializeParticles routine we set the positions of each particle by creating an instance of MdPointParticle and then adding it to the system. At this stage, it would be logical that we could set particle velocities as well.

We first create an instance of the random number generator class:

        StdLibErand48 *mRandomNumberGenerator = new StdLibErand48();

We can then use the the random number generator to create a vector of random numbers within the specific range we want, here within +/- sqrt(3.0):

        gvector v = mRandomNumberGenerator -> ComputeRandomVecT(-sqrt(3.0), sqrt(3.0));

Just as we used the SetPosition? routine, we can call the SetVelocity? routine to add the "v" vector.

        p -> SetVelocity(v);

This can be summarized int the following code exercpt:

void InitializeParticles()
{
        ipstream particle_pipe("mkcrystal -n=200 -d=0.85 -S");
        StdLibErand48 *mRandomNumberGenerator = new StdLibErand48();

        gvector x;         
        while(particle_pipe >> x)
        {
                MdPointParticle *p = new MdPointParticle;                  
                p -> SetPosition(x);
                gvector v = mRandomNumberGenerator -> ComputeRandomVecT(-sqrt(3.0), sqrt(3.0));
                p -> SetVelocity(v);
                AddParticle(p);
        }
        particle_pipe.close();
}

Completed Code

The following code creates particles on a simple cubic lattice using the "mkcrystal" routine, and additionally we use the random number generator class to create a random vector to set our velocities.

#include <glotzilla++.h>

class LennardJonesSimulationLattice : public MdSimulation
{
   public:
        LennardJonesSimulationLattice()
        {                                        
                InitializeParticles();
                SetInteraction(new LennardJones);
                SetForceRoutine(new BruteForce);
                SetBoundaryConditions(new PeriodicBoundary(10));
                SetIntegrationScheme(new VelocityVerletIntegrator);             
        }
        void InitializeParticles()
        {
                ipstream particle_pipe("mkcrystal -n=200 -d=0.85 -S");
                StdLibErand48 *mRandomNumberGenerator = new StdLibErand48();

                gvector x;         
                while(particle_pipe >> x)
                {
                        MdPointParticle *p = new MdPointParticle;                  
                        p -> SetPosition(x);
                        gvector v = mRandomNumberGenerator -> ComputeRandomVecT(-sqrt(3.0), sqrt(3.0));
                        p -> SetVelocity(v);
                        AddParticle(p);
                }
                particle_pipe.close();
        }
};

int main(int argc, char ** argv) 
{
        LennardJonesSimulationLattice my_simulation;
        opstream vis_pipe("vis3d");

        for(int timesteps=0; timesteps<10000; timesteps++)
        {
                my_simulation ++;
                if(timesteps%5 == 0)
                        my_simulation.PrintVisual(vis_pipe);
        }
        return 0;
}