# How can I write a unit test for this routine?

The first routine instantiates numAtoms number of atoms and distributes them in a 6*6*6 grid.

The second routine populates the atoms with initial velocities.

I want to write a unit test for this routine.

Can you give me some ideas on how to do it?

NOTE: I know how to write unit tests using GTest. That is not the issue. The issue is that I don't know what the expected values and actual values will be in this case, as the computations are pretty complicated.

void initializeCoordinates(std::vector<Atom>& atoms, int numAtoms)
{
const double a = SIGMA / 2.0; // Length unit a determined by σ = 2a

const double cellPositions[4][3] = {
{0.5 * a, 0.5 * a, 0.5 * a},
{0.5 * a, 1.5 * a, 1.5 * a},
{1.5 * a, 0.5 * a, 1.5 * a},
{1.5 * a, 1.5 * a, 0.5 * a}
};

const double cellShifts[3] = {2.0 * a, 2.0 * a, 2.0 * a};

int atomCount = 0;

for (int x = 0; x < numAtoms / 4; ++x) {
for (int y = 0; y < numAtoms / 4; ++y) {
for (int z = 0; z < numAtoms / 4; ++z) {
for (int i = 0; i < 4; ++i) {
double shiftX = x * cellShifts[0];
double shiftY = y * cellShifts[1];
double shiftZ = z * cellShifts[2];

if (atomCount >= numAtoms) {
// If the desired number of atoms is reached, return from the function
return;
}

Atom atom;
atom.number = atomCount + 1;
atom.coord_x = cellPositions[i][0] + shiftX;
atom.coord_y = cellPositions[i][1] + shiftY;
atom.coord_z = cellPositions[i][2] + shiftZ;

atoms.push_back(atom);

atomCount++;
}
}
}
}

// Move the system so that the center of mass is located at (0.0, 0.0, 0.0)
double centerOfMass[3] = {0.0, 0.0, 0.0};
double totalMass = 0.0;

for (const Atom& atom : atoms) {
centerOfMass[0] += atom.mass * atom.coord_x;
centerOfMass[1] += atom.mass * atom.coord_y;
centerOfMass[2] += atom.mass * atom.coord_z;
totalMass += atom.mass;
}

centerOfMass[0] /= totalMass;
centerOfMass[1] /= totalMass;
centerOfMass[2] /= totalMass;

for (Atom& atom : atoms) {
atom.coord_x -= centerOfMass[0];
atom.coord_y -= centerOfMass[1];
atom.coord_z -= centerOfMass[2];
}
}

void initializeVelocity(double T0, std::vector<Atom>& atoms)
{
#ifndef DEBUG
srand(time(NULL));
#endif
double totalMass = 0.0;
double vSquareSum = 0.0;

for (Atom& atom : atoms) {
atom.velocity_x = -6.0 + (rand() * 12.0) / RAND_MAX;
atom.velocity_y = -6.0 + (rand() * 12.0) / RAND_MAX;
atom.velocity_z = -6.0 + (rand() * 12.0) / RAND_MAX;
vSquareSum += atom.velocity_x * atom.velocity_x + atom.velocity_y * atom.velocity_y + atom.velocity_z * atom.velocity_z;
totalMass += atom.mass;
}

const double scaleFactor = sqrt((2.0 * T0 * K_B * totalMass) / vSquareSum);

for (Atom& atom : atoms) {
atom.velocity_x *= scaleFactor;
atom.velocity_y *= scaleFactor;
atom.velocity_z *= scaleFactor;
}

double centerOfMassVelocity[3] = {0.0, 0.0, 0.0};

for (const Atom& atom : atoms) {
centerOfMassVelocity[0] += atom.mass * atom.velocity_x;
centerOfMassVelocity[1] += atom.mass * atom.velocity_y;
centerOfMassVelocity[2] += atom.mass * atom.velocity_z;
}

centerOfMassVelocity[0] /= totalMass;
centerOfMassVelocity[1] /= totalMass;
centerOfMassVelocity[2] /= totalMass;

for (Atom& atom : atoms) {
atom.velocity_x -= centerOfMassVelocity[0];
atom.velocity_y -= centerOfMassVelocity[1];
atom.velocity_z -= centerOfMassVelocity[2];
}
}
$$$$
`

Divide and conquer

1. refactor your code, so that you have smaller methods with well defined purposes. Give these methods short descriptive names of what they are doing. You make the code more testable if these methods return numeric values, which you can easily interpret. Examples: Move your instantiation of the cell positions into getCellPositions(...) returning their positions, move the bulk of your iteration to a method (possibly marked with inline) "std::vector& atoms placeAtoms(...), move the calculation of the center of mass to a dedicated method "double[3] getCenterOfMass(..)".

2. for each of these methods, think of one regular input and the expected output you want. Write tests for these cases and structure them into three blocks: setup, call, assertions. Name these Tests with equally expressiveNames: "testCenterOfMass", "testKineticEnergyPlausible", "testNoOfAtoms","testAvgDistances",...

3. for each of the implementation methods, think of corner cases: null arguments, -infty, +infty, wrong array sizes etc. Does the method under test behave the way you want? Does it give a descriptive error in these corner cases? (i.e. invalidInput, overflow, underflow)

4. think about physical properties you expect to be met. These can be tested on a larger scope than individual methods. These could be the expected total mass, the expected kinetic energy of the system or the temperature. Minimal distance between two atoms (if two are in the same position you are in trouble:-)). When simulating a gass or a crystal you can find other physical properties in the literature that you want to reproduce like phase transitions.

5. Make use of the language features which help correctness. In c++ you can use the const keyword in various places to express wheter a method is allowed to alter an input parameter or not, or if a field is allowed to be changed at runtime. This helps the compiler with optimizations, but can also prevent bugs by being more expressive in your intent. Since C++17 you can even express certain calculations as constexpr indicating that they may already be performed by the compiler and will give compile-time-constant-results. Many of the methods after refactoring can be marked as constexpr, and you can make the compiler check it for you (similar to a unit test). The more you make use of these features, the more the compiler will help in 'testing' your code for correctness.

Other thoughts:

-consider posting this question to the dedicated code-review stack.

-Some random number generators produce identical output with identical seeds, you can make use of that fact for unit-testing without the need to alter your program.