Search code examples
pythonnumpyphysicsvpython

Sodium chloride crystal


I'm trying to do this exercise:

A sodium chloride crystal has sodium and chlorine atoms arranged on a cubic lattice but the atoms alternate between sodium and chlorine, so that each sodium is surrounded by six chlorines and each chlorine is surrounded by six sodiums. Create a visualization of the sodium chloride lattice using two different colors to represent the two types of atoms.

My code is:

from vpython import *
from numpy import *

L = 5
R = 0.5
for i in range(-L,L,2):
    for j in range(-L,L,2):
        for k in range(-L,L,2):
            sphere(pos=vector(i,j,k),radius = R, color = vec(0,1,1))

for l in range(-L+1,L+1,2):
    for m in range(-L+1,L+1,2):
        for n in range(-L+1,L+1,2):
            sphere(pos=vector(l,m,n),radius = R, color = vec(1,1,0))

But I get this figure: enter image description here

From this angle, you can see that it is not correct because there are columns of the same color. What I am doing wrong?


Solution

  • First off, the provided code does not generate a well aligned crystal lattice for Sodium Chloride with the expected orientation.

    Consider the first few numbers of each of the output.

    >>> for i in range(-L,L,2):
    ...     for j in range(-L,L,2):
    ...         for k in range(-L,L,2):size = 5
    ...             print(i, j, k)
    ... 
    -5 -5 -5
    -5 -5 -3
    -5 -5 -1
    ...
    >>> for l in range(-L+1,L+1,2):
    ...     for m in range(-L+1,L+1,2):
    ...         for n in range(-L+1,L+1,2):
    ...             print(l, m, n)
    ... 
    -4 -4 -4
    -4 -4 -2
    -4 -4 0
    ...
    

    If the former spheres offsets the latter spheres by one unit, every sphere would have instead 10 neighbors instead of 6. This explains why the spheres are rendered rather far apart and that they appear to have 10 neighbors each, thus what was attempted does not appear as a tightly packed cubic lattice but rather all spread out. Effectively, what you did wrong was applying an offset to the lattice when you really shouldn't, introducing the gaps between the spheres to make this look visually off.

    In order to do this correctly, take this slowly and one step at a time, and iterate every cell at that layer and render the sphere in one go. Let's start with a 2-D checkerboard pattern first.

    Consider the following instead:

    from vpython import *
    from numpy import *
    
    size = 5
    radius = 0.5
    # predefine the colours
    elements = [
        vec(0, 1, 1),
        vec(1, 1, 0),
    ]
    
    z = 0  # pin z to just a single layer 
    for y in range(-size, size):
        for x in range(-size, size):
            sphere(pos=vector(x, y, z), radius=radius, color=elements[(x + y) % 2])
    

    The color is calculated by simply adding the x and y, then take the modulus 2 of the result to determine which color to use, as sums of an odd with even number will produce an odd number (otherwise they are even), and this mathematical property when rendered out will have our first layer. A related thread on this here on StackOverflow.

    Turns out this property also applies with the Z layer

    for z in range(-size, size):
        for y in range(-size, size):
            for x in range(-size, size):
                sphere(pos=vector(x, y, z), radius=radius, color=elements[(x + y + z) % 2])
    

    That will produce what you expected. Now if you want to modify the radius argument and color together, you might want to specify the elements as two dictionaries that can be passed in via keyword argument expansion:

    size = 5
    elements = [
        {'radius': 0.6, 'color': vec(0, 1, 0)},
        {'radius': 0.4, 'color': vec(0.6, 0.1, 1)},
    ]
    
    for z in range(-size, size):
        for y in range(-size, size):
            for x in range(-size, size):
                sphere(pos=vector(x, y, z), **elements[(x + y + z) % 2])                                                                             
    

    That would mimic the crystal structure diagram for Sodium Chloride that might be found on Wikipedia (view rotated to show the result more clearly).

    Generated diagram