Note: This is a follow-up to my post on the Madelung constant (it’s for NaCl). If you don’t know what that is yet, take a look at it here.
I wrote the Madelung constant finder code in Python, so you’ll need some interpreter if you want to run it alongside this explanation.
If you don’t have it already, I recommend downloading the Anaconda distribution for Python. It comes preloaded with a bunch of useful packages (NumPy, SciPy, Matplotlib, etc.) and Jupyter, an IDE good for demonstrating what your code can do in an organized and runnable environment. I’m uploading the Jupyter notebook (159 downloads ) of this code as well as the raw .py scripts (176 downloads ) – I highly recommend the Jupyter option if you want to follow along closely.
Coordinate Generator
The first thing I did was establish the best logical framework for computing the constant. In my opinion, the best option is to generate a list of lists—with each nested list representing a “triplet” coordinate set—and then calculate the influence for each (Not like real triplets, since the things inside aren’t identical. I just call them that in my syntax, so I’ll use the name here too for consistency). You might be quick to point out the potential for some shorter, more elegant, and vectorised solutions here, but the triplet list method is suprisingly fast (especially if you only generate a few and multiply, as I demonstrate) and it carries the benefit of being really easy to understand just from looking at the code.
So, how do you generate a list of all possible 3-d coordinate sets? You don’t. Well, you can, but it’s slower. You’ll see what I’m talking about soon.
My friend whipped up a script for generating coordinate sets, which I refined to stop at certain shells (you just end it when the maximum number in a set is past your desired shell number). It looks like this:
while (triplets[-1][0 and 1 and 2]) < shell_number:
lasttrip = (triplets[-1])
if lasttrip[0] == lasttrip[1]:
if lasttrip[1] == lasttrip[2]:
triplets.append([0, 0, lasttrip[2]+1])
elif lasttrip[1] < lasttrip[2]:
triplets.append([0, lasttrip[1]+1, lasttrip[2]])
elif lasttrip[0] < lasttrip[1] and lasttrip[1] <= lasttrip[2]:
triplets.append([lasttrip[0]+1, lasttrip[1], lasttrip[2]])
The first line makes it loop through the rest until any of the coordinates in the last generated set are more than the shell number. From there, it goes through a bunch of conditional statements that determine what the next coordinate set should be based on the last generated one—sometimes adding one to a term, other times flipping a term to zero and adding one, or rarely, flipping the whole thing (like from [1,1,1] to [0,0,2]).
You may notice something strange about this code. Namely, that it doesn’t actually generate every possible coordinate set. If you run this with, say, shell_number = 3, you’ll get this as the output into triplets (to run this snippet yourself, you’ll first need to set triplets as a list of lists, with its first term as [0,0,1]):
[0, 0, 1] [0, 1, 1] [1, 1, 1] [0, 0, 2] [0, 1, 2] [1, 1, 2] [0, 2, 2] [1, 2, 2] [2, 2, 2] [0, 0, 3] [0, 1, 3] [1, 1, 3] [0, 2, 3] [1, 2, 3] [2, 2, 3] [0, 3, 3] [1, 3, 3] [2, 3, 3] [3, 3, 3] ...
They sort of just count up, as if it were some kind of weird binary on the first shell, trinary on the second, and so on. This, of course, misses all the negative numbers and position combinations of numbers in between. But we can solve this by multiplying the influence of each triplet by the amount of atoms that would be at the same distance. Since direction doesn’t matter (only distance and amount), it’s effectively the same thing.
Equal Distance Set Generator
The logic tree for doing this is pretty long, but it’s also pretty interesting. In code form, it looks like this:
for i in range(len(triplets)):
coordset=triplets[i]
if coordset[0] == 0:
if coordset[1] == 0:
eqdist = 6
elif coordset[1] == coordset[2]:
eqdist = 12
else:
eqdist = 24
elif coordset[0] == coordset[1]:
if coordset[1] == coordset[2]:
eqdist = 8
else:
eqdist = 24
else:
if coordset[1] == coordset[2]:
eqdist = 24
else:
eqdist = 48
eqdistset.append(eqdist)
This makes a list the same length as triplets that’s filled with the number of atoms equidistant to its corresponding triplets entry.
This list is also the number of atoms on each sphere of increasing radii (I touched on this sequence a bit in my original post. It’s also on the online encyclopedia of integer sequences). It was pretty fun to work out how many atoms were at the same distance based on some general facts about the particular coordinate set. If you have time, I would highly recommend trying it out yourself for fun. It’ll also give you more insight into the weirdness of that sequence.
Denominator Generator
I also wrote in my last post that, in order to converge faster, we needed to account for the fraction of the atom in the shell. To picture how this works, imagine each shell as the surface of a cube with its corners centered on an atom. The atoms at those corners only have 1/8 on the inside of the cube, the edge atoms have 1/4, and face atoms 1/2.
To write code for this, I first thought you would need to add the fraction of each layer and then store the remainder to add onto the next shell up. But if you really think about, you only need to fractionalize the last layer. For everything in between, it’s as if you’ve already added both the inside and outside parts when you just don’t do anything to them.
This code generates a list of denominators for the fraction in the shell of the corresponding atom in triplets:
for i in range(len(triplets)):
coordset = triplets[i]
if coordset[0] == coordset[1]:
if coordset[1] == coordset[2]:
denom = 8
else:
denom = 2
elif coordset[1] == coordset[2]:
denom = 4
else:
denom = 2
denomlist.append(denom)
Short and sweet. It’s just checking if a set represents a corner, edge, or face atom and then appending a value to denomlist accordingly.
Putting it All Together
Now let’s do the sum. We’ll generate a list of addends based on the other lists we made. The generating algorithm needs to determine if the coordinate set represents a sodium or chlorine atom so we know whether to add or subtract (a simple way to do this is to add the coordinate values together—if it’s odd then the atom is opposite the reference, even being the same). It then needs to determine if the coordinate set represents an atom on the edge of the shell, so it knows whether or not to use denomlist.
Finding the distances is a trivial task with numpy:
triplets = np.array(triplets) distset = np.sqrt(np.sum(triplets**2,axis = 1))
And so the final generator looks like this:
for i in range(len(triplets)):
coordset=triplets[i]
if abs(coordset[0]+coordset[1]+coordset[2])%2 == 1:
if max(coordset) == shell_number:
addlist.append(eqdistset[i]/(denomlist[i]*distset[i]))
else:
addlist.append(eqdistset[i]/distset[i])
else:
if max(coordset) == shell_number:
addlist.append(-eqdistset[i]/(denomlist[i]*distset[i]))
else:
addlist.append(-eqdistset[i]/distset[i])
Now we just need to extract our data.
print("Final calculated Madelung constant for NaCl is:", sum(addlist))
And we’re done!