ParM  parm
A molecular dynamics library


Based on LJatoms.cpp


Include the necessary header files:

Make an AtomVec

Create 20 atoms, with mass 1.0:

uint Natoms = 20;
boost::shared_ptr<AtomVec> atomptr(new AtomVec(Natoms, 1.0));

We use boost::shared_ptr for the major objects, to prevent memory leaks / segmentation faults and integrate more smoothly with Python. If you are not familiar with shared_ptr, its a wrapper around a pointer such that multiple copies can be made, and the object will be deleted when the last shared_ptr to it is deleted. See a reference for more details.

Make a Box

const flt L = 20;
boost::shared_ptr<OriginBox> obox(new OriginBox(L));

Here, we'll use a periodic box with sides of length L.

Make at least one Interaction

We'll use the Lennard-Jones Interaction, \(V\left(r\right)=\varepsilon\left(1-\frac{\sigma^{6}}{r^{6}}\right)^{2}\), where \(\sigma\) is the size of each particle.

For computational purposes, we "truncate and shift" this potential at \(2.5\sigma\), a standard practice. We will also use a Neighbor list, which keeps track of which atoms are in the vicinity of each other.

We create that Interaction, and give its NeighborList a variable name:

const flt sigma = 1.;
const flt sigcut = 2.5;
boost::shared_ptr<NListed<EpsSigCutAtom, LennardJonesCutPair> >
LJ(new NListed<EpsSigCutAtom, LennardJonesCutPair>(obox, atomptr, 0.4*(sigcut*sigma)));
boost::shared_ptr<NeighborList> neighbor_list = LJ->neighbor_list();
Note the use of templating: `NListed<A, P>` is a template for an `Interaction` that *can* be
used with a NeighborList, and can be used with a number of different potentials, or custom ones.
#### Initialize the Atoms
for (uint i=0; i < atomptr->size(); i++){
atoms[i].x = obox->rand_loc(); // random location in the box
atoms[i].v = rand_vec(); // from a Gaussian
atoms[i].f = Vec();
atoms[i].a = Vec();
// Add it to the Lennard-Jones potential
LJ->add(EpsSigCutAtom(epsilon, sigma, atoms.get_id(i), sigcut));
// Update the neighbor list

Now we have given every Atom a position and a random velocity.

Make a Collection

A Collection is a class that represents an integrator. There are various kinds of integrators: constant energy, constant temperature, energy minimization, etc. Here, we choose a Velocity Verlet integrator, which is constant energy:

const flt dt = 0.0001;
CollectionVerlet collec = CollectionVerlet(boost::static_pointer_cast<Box>(obox), atomptr, dt);

Link the Interaction, NeighborList, and Collection

We made our Collection, but it needs to know what potential to use with the atoms, and that it needs to update the NeighborList:


Running the simulation

Let's run it for 500,000 timesteps, outputting the current energy, kinetic_energy energy, and potential energy as we go:

for(uint i=0; i<500; i++){
for(uint j=0; j<1000; j++){
writefile(outfile, atoms, *obox);
cout << (500-i) << " E: " << << " K: " << collec.kinetic_energy()
<< " U: " << LJ->energy(*obox) << "\n";

That's all! To see this example fully fleshed out, see LJatoms.cpp