Calculating x,y,z positions in 3D space to create face-centered cubic lattice I am a university student and am trying to find a way to calculate the (x,y,z) positions for an arbitrary number coordinates along a face centered cubic lattice, centered around the origin. 
I am looking for an equation to calculate the positions (x, y, z) of n number of points that "wraps" around the origin. For example, if I want to plot 12 points, each should be sqrt(2) away from the origin and the corresponding points should be (0,1,1),(1,0,1),(1,1,0), (0,−1,1),(−1,0,1),(−1,1,0), (0,1,−1),(1,0,−1),(1,−1,0), (0,−1,−1),(−1,0,−1),(−1,−1,0), but what if I want to calculate for 100 points? 500 points? I am trying to create an equation that does just that.
I know I need to utilize the properties of a face-centered cubic lattice but am lost on how to approach this problem. I was thinking about using the corresponding Gram Matrix to aid transformations but do not know if this is the correct approach. 
Any answers/advice/guidance would be greatly appreciated. Thank you.
 A: Hint: Integer points with $x+y+z$ are even.
A: So, you want to generate the 3D coordinates for $N$ lattice points, for the $N$ lattice points within distance $r$ from origin, in a face-centered cubic lattice, in order of distance from origin. 
If we look at the lattice cell at $(x_0, y_0, z_0)$ (i.e., $x_0 \le x \lt x_0 + 1$, $y_0 \le y \lt y_0 + 1$, and $z_0 \le z \lt z_0 + 1$), we find that there are lattice sites at coordinates
$$\begin{array}{lll|l}
x & y & z & \text{Description} \\
\hline
x_0     & y_0     & z_0     & \text{Corner} \\
x_0 + \frac{1}{2} & y_0 + \frac{1}{2} & z_0               & \text{Face} \\
x_0 + \frac{1}{2} & y_0               & z_0 + \frac{1}{2} & \text{Face} \\
x_0               & y_0 + \frac{1}{2} & z_0 + \frac{1}{2} & \text{Face} \\
\end{array}$$
The distance from origin to each such lattice point is
$$\begin{array}{rll}
d_0 = & \; & \sqrt{x_0^2 + y_0^2 + z_0^2} \\
d_1 = & \sqrt{\left(x_0 + \frac{1}{2}\right)^2 + \left(y_0 + \frac{1}{2}\right)^2 + z_0^2} = & \sqrt{x_0^2 + y_0^2 + z_0^2 + x_0 + y_0 + \frac{1}{2}} \\
d_2 = & \sqrt{\left(x_0 + \frac{1}{2}\right)^2 + y_0^2 + \left(z_0 + \frac{1}{2}\right)^2} = & \sqrt{x_0^2 + y_0^2 + z_0^2 + x_0 + z_0 + \frac{1}{2}} \\
d_3 = & \sqrt{x_0^2 + \left(y_0 + \frac{1}{2}\right)^2 + \left(z_0 + \frac{1}{2}\right)^2} = & \sqrt{x_0^2 + y_0^2 + z_0^2 + y_0 + z_0 + \frac{1}{2}} \\
\end{array}$$
Because of this, the distance to each lattice point $i$ is $\sqrt{k_i + \frac{1}{2}}$, where $0 \le k_i \in \mathbb{Z}$ for all $i$.
OP is using a coordinate system where all sites have integer coordinates: "corner" sites occur at even $x$, $y$, and $z$ coordinates, and "face" sites have two odd coordinates and one even coordinate. Thus, $x + y + z$ is even for all lattice points, and only for lattice points. In this case, the distance to each lattice point is $\sqrt{2 k_i}$, where $0 \le k_i \in \mathbb{Z}$ for all $i$.
The number of points with integer coordinates at distance $\sqrt{d}$ from origin, is an open problem. In 2D, it is known as Gauss circle problem. In this case, we have a related problem one could call "Gauss sphere problem using a face-centered cubic lattice".
OP wants to construct an equation, most likely a vector-valued function $\vec{f}(i)$, or possibly an integer vector sequence $\vec{f}_i$, where $i$ is an integer index, and $\vec{f}$ yields the location of the $i$'th lattice site, when the lattice sites are sorted according to their distance to origin.  Even if we assume that the order of sites having the same distance is irrelevant, this is an open problem: No known solutions exist. 
If we use OP's coordinate system, where all lattice points have integer coordinates ($x + y + z$ is even for all lattice sites), then the sequence $n_d$ describing the number of FCC lattice sites at distance $\sqrt{2d}$ from origin is OEIS A004015; $0 \le d \in \mathbb{Z}$. The sequence $N_d$ describing the number of FCC lattice sites within distance $\sqrt{2d}$ from origin, $N_d = \sum_{i=0}^d n_i$, is OEIS A119869.
(That also means that the number of lattice sites $N$ OP is interested in is basically $N_d$, with $N_{d-1} \lt N \le N_d$, assuming $N_i = -1 \forall i \lt 0$.)

What we can do, however, is construct numerical algorithms for generating the lattice site coordinates for sites at or within distance $\sqrt{d}$ from origin; their main difference being the amount of memory used. The naïve approach has something like $O(d^3)$ time complexity, but it can be somewhat reduced by clever approaches or by using memory to store details of the FCC lattice.
In a computer program where I needed each lattice site set $f_{d, i}$ (lattice sites at distance $\sqrt{d}$ from origin) separately, I would probably use a purely numeric approach with a fixed memory use (except for the lattice site coordinate sets).
The idea is that lattice sites at distance $\sqrt{d}$ from origin have $$\left\lfloor \sqrt{d - x^2 - y^2} \right\rfloor \le z \le \left\lceil \sqrt{d - x^2 - y^2} \right\rceil$$
so only a small subvolume of the positive octant needs to be examined to locate the lattice sites.
This seems to be approximately $O(d)$, i.e. linear in the distance squared; or $O(r^2)$ where $r$ is the distance.
Here is a very simple Python 3 implementation:

from math import sqrt, ceil, floor

def fcc_sites(dsqr):
    """Generates integer FCC lattice sites at distance sqrt(dsqr) from origin.
       Use even dsqr.""""

    # There are no sites at odd dsqr.
    if int(dsqr) % 2 == 1:
        return

    xmax = int(ceil(sqrt(dsqr)))

    # x loops from 0 to xmax, inclusive:
    for x in range(xmax+1):
        try:
            ymax = int(ceil(sqrt(dsqr - x*x)))
        except:
            ymax = -1

        for y in range(ymax + 1):
            try:
                ztmp = sqrt(dsqr - x*x - y*y)
            except:
                ztmp = 0

            zmax = int(ceil(ztmp))
            zmin = int(floor(ztmp))
            zz = dsqr - x*x - y*y
            for z in range(zmin, zmax+1):
                if zz == z*z and (x + y + z) % 2 == 0:
                    if x == 0 and y == 0 and z == 0:
                        yield((0, 0, 0))
                    elif x == 0 and y == 0:
                        yield((0, 0, z))
                        yield((0, 0, -z))
                    elif x == 0 and z == 0:
                        yield((0, y, 0))
                        yield((0, -y, 0))
                    elif y == 0 and z == 0:
                        yield((x, 0, 0))
                        yield((-x, 0, 0))
                    elif x == 0:
                        yield((0, y, z))
                        yield((0, -y, z))
                        yield((0, y, -z))
                        yield((0, -y, -z))
                    elif y == 0:
                        yield((x, 0, z))
                        yield((-x, 0, z))
                        yield((x, 0, -z))
                        yield((-x, 0, -z))
                    elif z == 0:
                        yield((x, y, 0))
                        yield((-x, y, 0))
                        yield((x, -y, 0))
                        yield((-x, -y, 0))
                    else:
                        yield((x, y, z))
                        yield((-x, y, z))
                        yield((x, -y, z))
                        yield((-x, -y, z))
                        yield((x, y, -z))
                        yield((-x, y, -z))
                        yield((x, -y, -z))
                        yield((-x, -y, -z))

One could make this even more efficient, by utilizing the symmetries even further. For example, one could consider only sites with $0 \le x \le y \le z$, and permute the coordinates in addition to mirroring. Or, one could construct a walker, essentially one eighth of a circle in the $yx$ plane, again using symmetries; based on midpoint circle or Bresenham circle algorithms.
Also note that to generate the lattice points within some distance from origin, without sorting them according to distance, can be done much more efficiently, as only the "surface cells'" lattice sites need to be examined for their distance, and all the inner cells can simply be generated without examination.
The above Python code is only lightly tested, so if you find any issues with it, please let me know in a comment so I can verify and fix.
A: I'd like to just point you to the topic of Waterman polyhedra. He'd got the idea to consider polyhedra as being the convex hull of those lattice points inside some given radius around some given lattice point.
The wiki page does consider only the CCP: https://en.wikipedia.org/wiki/Waterman_polyhedron
Here is his own home page, showing up different lattices: http://watermanpolyhedron.com/polyhedra.html
--- rk
A: For the sake of completeness, I wanted to follow up with the solution to my problem. For clarity, my question was more precisely: "how to translate lattice coordinates to Cartesian coordinates." 
I have defined the labeling scheme for a face-centered cubic lattice according to the following picture. Each vertex in the lattice has twelve neighbors and this is how I locally defined the lattice 'directions' from a vertex centered at (0,0,0). 

This is my thought process for solving this problem:
Within the same z-plane, (for example, neighbors 0-5 are in the same z-plane) the distance between the neighboring vertices is equal to 1 & any three neighboring particles in the same z-plane form an equilateral triangle with a tetrahedral hole. Between adjacent z-planes, the stack of two sets of three vertices creates an octahedral hole.  (https://chem.libretexts.org/Bookshelves/Inorganic_Chemistry/Supplemental_Modules_(Inorganic_Chemistry)/Crystal_Lattices/Solids/Tetrahedral_and_Octahedral_Sites). 
From the picture, vertex #6 is located at (0,0,1) but what is the coordinate in Cartesian space? vertex #6 is located directly above the center of the equilateral triangle that is created by vertex located at (0,0,0), vertex #0 and vertex #1. After a few calculations, vertex #6's cartesian coordinate is (.5, -sqrt(3.0)/4.0, Z). The Z-component of the Cartesian coordinate now is the only unknown... Since vertex #6 sits in the middle of an equilateral triangle, the four particles create a tetrahedron and the z-component is the height of that tetrahedron. After some geometric calculations, I have come up with the following equation to translate from lattice to Cartesian coordinates.
double x = i + j/2.0 + k/2.0;
double y =  sqrt(3.0)/2.0 * j - sqrt(3.0)/4.0 * k;
double z =  sqrt(6.0)/3.0 * k;

