putting particles on fcc lattice

[PhysicsIsFun]

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 a link to the picture of the lattice structure:
https://www.physics-in-a-nutshell.com/img/content/solid-state-physics/fcc-coordination-number-nearest-neighbours.svg

The lattice is spanned by the three lattice vectors
a1= a/2 (0,1,1), a2=a/2 (1,0,1), a3=a/2 (1,1,0)

How would you do this most easily?

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.

edit: mb I do it layer wise and forget about forming the lattice vector combinations all together...

Regards

[PhysicsIsFun]

I think I have my solution.
Just in case someone is curious^^

#include <fstream>
#include <string>
#include <vector>
//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;

vector<double> x0_values(columns), y0_values(rows), z0_values(layers),
               x1_values(columns), y1_values(rows), z1_values(layers);

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

        z0_values.at(i)=z0;
        z1_values.at(i)=z1;

        x0=-a;

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

            y0=-a;

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

                y0_values.at(k)=y0;
                y1_values.at(k)=y1;

                fcc << x0_values.at(j) << "\t" << y0_values.at(k) << "\t" << z0_values.at(i) << endl;
                fcc << x1_values.at(j) << "\t" << y1_values.at(k) << "\t" << z0_values.at(i) << endl;
                fcc << x0_values.at(j) << "\t" << y1_values.at(k) << "\t" << z1_values.at(i) << endl;
                fcc << x1_values.at(j) << "\t" << y0_values.at(k) << "\t" << z1_values.at(i) << endl;
            }
        }
}

fcc.close();
}

edit: vectors are not even necessary here, I guess. So better version:

#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();
}