Setting up an FCC lattice in a simulation

[SchroedingersLion]

Greetings,

for a Monte Carlo simulation of N hard spheres, I want to set up the initial configuration as a fcc lattice. I am thinking about how to do it most efficiently.

Here is the lattice structure:

The 3 lattice vectors, each one pointing to a different type of the 12 nearest neighbors (see colors in picture), are given by
$$\vec a_{1}=\frac a 2 \begin{pmatrix} 0 \\ 1 \\ 1 \end{pmatrix},
\vec a_{2}=\frac a 2 \begin {pmatrix} 1 \\ 0 \\ 1 \end {pmatrix},
\vec a_{3}=\frac a 2 \begin {pmatrix} 1 \\ 1 \\ 0 \end {pmatrix}
$$

How would you do this most easily in C++?

My first idea was to simply set up the nearest neighbors. First, I use the three lattice vectors separately to set up the first 3 nearest neighbors. Then I combine two of them to get the other 9. Then I combine all 3 to get the next-nearest neighbors.

But what then? I have to repeat these steps with different pre-factors in the linear combinations, but I don't know how to do it systematically in a loop.Regards,
SL

 

[SchroedingersLion]

Ok, in case someone needs this as well, here is my routine:

#include <fstream>
#include <string>

//function that writes the coordinates of fcc lattice sites with given
//lattice constant, numbers of rows, columns and layers to an output file

using namespace std;

void generate_fcc(double a, int rows, int columns, int layers){

ofstream fcc("fcc_coordinates_a=" + to_string(a) + "_" +
             to_string(rows) + "rows_" + to_string(columns) + "columns_" + to_string(layers) + "layers.dat");

//necessary initializations to get coordinates of 1st particle after incrementations at beginning of loops
double x0=-a, y0=-a, z0=-a, x1, y1, z1=-a/2;

// print coordinates, layer by layer, columns by columns, row by row
for (int i=0; i<layers; i++){
        z0+=a;
        z1+=a;

        x0=-a;

        for (int j=0; j<columns; j++){
            x0+=a;
            x1=x0 + a/2;

            y0=-a;

            for (int k=0; k<rows; k++){
                y0+=a;
                y1=y0 + a/2;

                fcc << x0 << "\t" << y0 << "\t" << z0 << endl;
                fcc << x1 << "\t" << y1 << "\t" << z0 << endl;
                fcc << x0 << "\t" << y1 << "\t" << z1 << endl;
                fcc << x1 << "\t" << y0 << "\t" << z1 << endl;
            }
        }
}

fcc.close();
}