I’m currently in the process of porting over my old WordPress blog to GitHub Pages, which supports Jekyll. I was getting tired of the bloat dealing with Wordpress, and wanted a cleaner, snappier website.
That and now I can blog in markdown, which is a great advantage to using Jekyll! From the terminal! And version control and hosting with git and GitHub!
I’m a little late to the party, but glad I made it.
I released some new code today that should allow you to easily rank NFL teams for any given season (back through 2009) using either the Colley or Massey methods. Both methods have several weighting schemes implemented, and with the Massey method you can additionally rank according to point differential or total yard differential.
You can check it out on my GitHub here.
I’ve become increasingly dissatisfied with my old NFL Colley ranking system and its derivatives, mostly because (a) I didn’t write it in an object-oriented style, and (b) it became very (VERY) hard to modify later on. The new structure, by encapsulating everything in a “Season” object should keep everything pretty self-contained.
For example, now
from season import Season
season = Season()
season.year = 2014 # default season is 2015
season.massey() # Massey ranking method, unweighted
for team in season.rating:
print team Which gives the output
['NE', 11.407605241747872]
['DEN', 9.6038904227178108]
['SEA', 9.5169656599013628]
['GB', 8.2526935620363258]
...
['OAK', -8.9899785292554917]
['TB', -9.8107991356959232]
['JAC', -10.427123019821691]
['TEN', -11.805248019821692]So obviously the NE Patriots were ranked #1 for the 2014 season with this method. You’ll recall they ended up winning the Super Bowl that same season.
So anyway, I’m starting over and making use of some great NFL APIs that I have found elsewhere on GitHub. In particular, I am using nflgame, which does a lot of the heavy lifting for me associated with scraping necessary data.
Check it out if this sounds like it may be something you’re interested in!
Recently, I found myself executing the same commands (or some variation thereof) at the command line.
Over and over and over and over.
Sometimes I just write a do loop (bash) or for loop (python). But the command I was executing was taking almost a minute to finish. Not terrible if you are doing this once, but several hundred (thousand?) times more and I get antsy.
If you are curious, I was generating .cube files for a large quantum dot, of which I needed the molecular orbital density data to analyze. There was no way I was waiting half a day for these files to generate.
So instead, I decided to parallelize the for loop that was executing my commands. It was easier than I thought, so I am writing it here not only so I don’t forget how, but also because I’m sure there are others out there like me who (a) aren’t experts at writing parallel code, and (b) are lazy.
Most of the following came from following along here.
First, the package I used was the joblib package in python. I’ll assume you have it installed, if not, you can use pip or something like that to get it on your system. You want to import Parallel and delayed.
So start off your code with
from joblib import Parallel, delayed If you want to execute a system command, you’ll also need the call function from the subprocess package. So you have
from joblib import Parallel, delayed
from subprocess import call Once you have these imported, you have to structure your code (according to the joblib people) like so:
import ....
def function1(...):
...
def function2(...):
...
...
if __name__ == '__main__':
# do stuff with imports and functions defined about
... So do you imports first (duh), then define the functions you want to do (in my case, execute a command on the command line), and then finally call that function in the main block.
I learn by example, so I’ll show you how I pieced the rest of it together.
Now, the command I was executing was the Gaussian “ cubegen” utility. So an example command looks like
cubegen 0 MO=50 qd.fchk 50.cube 120 h Which makes a .cube file (50.cube) containing the volumetric data of molecular orbital 50 (MO=50) from the formatted checkpoint file (qd.fchk). I wanted 120 points per side, and I wanted headers printed (120 h).
Honestly, the command doesn’t matter. If you want to parallelize
ls -lh over a for loop, you certainly could. That’s not my business.
What does matter is that we can execute these commands from a python script using the call function that we imported from the subroutine package.
So we replace our functions with the system calls
from joblib import Parallel, delayed
from subprocess import call
def makeCube(cube,npts):
call(["cubegen","1","MO="+str(cube),"qd.fchk",str(cube)+".cube",
str(npts),"h"])
def listDirectory(i): #kidding, sorta.
call(["ls", "-lh"])
if __name__ == '__main__':
# do stuff with imports and functions defined about
... Now that we have the command(s) defined, we need to piece it together in the main block.
In the case of the makeCube function, I want to feed it a list of molecular orbital (MO) numbers and let that define my for loop. So let’s start at MO #1 and go to, say, MO #500. This will define our inputs. I also want the cube resolution (npts) as a variable (well, parameter really).
I’ll also use 8 processors, so I’ll define a variable num_cores and set it to 8. Your mileage may vary. Parallel() is smart enough to handle fairly dumb inputs.
(Also, if you do decide to use cubegen, like I did, please make sure you have enough space on disk.)
Putting this in, our code looks like
from joblib import Parallel, delayed
from subprocess import call
def makeCube(cube,npts):
call(["cubegen","1","MO="+str(cube),"qd.fchk",str(cube)+".cube",
str(npts),"h"])
def listDirectory(i): #kidding, sorta
call(["ls", "-lh"])
if __name__ == '__main__':
start = 1
end = 501 # python's range ends at N-1
inputs = range(start,end)
npts = 120
num_cores = 8Great. Almost done.
Now we need to call this function from within Parallel() from joblib.
results = Parallel(n_jobs=num_cores)(delayed(makeCube)(i,npts)
for i in inputs) The Parallel function (object?) first takes the number of cores as an input. You could easily hard code this if you want, or let Python’s multiprocessing package determine the number of CPUs available to you. Next we call the function using the delayed() function. This is “a trick to create a tuple (function, args, kwargs) with a function-call syntax”.
It’s on the developer’s web page. I can’t make this stuff up.
Then we feed it the list defined by our start and end values.
If you wanted to list the contents of your directory 500 times and over 8 cores, it would look like (assuming you defined the function and inputs above)
results = Parallel(n_jobs=8)(delayed(listDirectory)(i)
for i in inputs) Essentially we are making the equivalence that
delayed(listDirectory)(i) for i in inputsis the same as
for i in inputs:
listDirectory(i)Does that make sense? It’s just
delayed(function)(arguments)instead of
function(arguments)Okay. Enough already. Putting it all together we have:
from joblib import Parallel, delayed
from subprocess import call
def makeCube(cube,npts):
call(["cubegen","1","MO="+str(cube),"qd.fchk",str(cube)+".cube",
str(npts),"h"])
def listDirectory(i): #kidding, sorta
call(["ls", "-lh"])
if __name__ == '__main__':
start = 1
end = 501 # python's range ends at N-1
inputs = range(start,end)
npts = 120
num_cores = 8
results = Parallel(n_jobs=num_cores)(delayed(makeCube)(i,npts)
for i in inputs) There you have it!
Lately I have been diagonalizing some nasty matrices.
Large. Non-Hermitian. Complex. Matrices. The only thing I suppose I have going for me is that they are relatively sparse.
Usually I haven’t have much of a problem getting eigenvalues. Eigenvalues are easy. Plug into ZGEEV, compute, move on.
The problem I ran into came when I wanted to use the eigenvectors. If you are used to using symmetric matrices all the time, you might not realize that non-Hermitian matrices have two sets of eigenvectors, left and right. In general, these eigenvectors of a non-Hermitian matrix are not orthonormal. But you can biorthogonalize them.
Why do you care if your eigenvectors are biorthogonalized?
Well, if your matrix corresponds to a Hamiltonian, and if you want to compute wave function properties, then you need a biorthonormal set of eigenvectors. This happens in linear response coupled cluster theory, for example. It is essential for a unique and physical description of molecular properties, e.g. transition dipole moments.
Now, with Hermitian matrices, your left and right eigenvectors are just conjugate transposes of each other, so it’s super easy to orthogonalize a set of eigenvectors. You can compute the QR decomposition (a la Gram-Schmidt) of your right eigenvectors \(\mathbf{C}\) to get
\[\begin{equation} \mathbf{C} = \mathbf{QR} \end{equation}\]where \(\mathbf{Q}\) is your set of orthogonalized eigenvectors. (Usually they are orthogonalized anyway during your eigenvalue decomposition.)
For non-Hermitian matrices, this situation is different. You have two eigenvalue equations you want to solve:
\(\mathbf{HR} = \mathbf{RE}\), and, \(\mathbf{LH} = \mathbf{LE}\)
where \(\mathbf{R}\) and \(\mathbf{L}\) are your right and left eigenvectors, and \(\mathbf{H}\) and \(\mathbf{E}\) are your matrix and eigenvalues. The eigenvalues are the same regardless of which side you solve for.
To biorthonormalize \(\mathbf{L}\) and \(\mathbf{R}\), we want to enforce the constraint
\[\mathbf{LR} = \mathbf{I}\]which says that the inner product of these two set of vectors is the identity. This is what biorthonormalization means.
Many times \(\mathbf{LR} = \mathbf{D}\) where \(\mathbf{D}\) is diagonal. This is easy. It’s already orthogonal, so you can just scale by the norm.
If that isn’t the case, how do you do it? One way would be to modify Gram-Schmidt. You could do that, but you won’t find a LAPACK routine to do it for you (that I know of). So you’d have to write one yourself, and doing that well is time-consuming and may be buggy. Furthermore, I’ve found the modified Gram-Schmdt to run into to serious problems when you encounter degenerate eigenvalues. In that case, the eigenvectors for each degenerate eigenvalue aren’t unique, even after constraining for biorthonormality, and so it’s tough to enforce biorthonormality overall.
Here’s a trick if you just want to get those dang eigenvectors biorthonormalized and be on your way. The trick lies in the LU decomposition.
Consider the following. Take the inner product of \(\mathbf{L}\) and \(\mathbf{R}\) to get the matrix \(\mathbf{M}\).
\[\mathbf{LR} = \mathbf{M}\]Now take the LU decomposition of \(\mathbf{M}\)
\[\mathbf{M} = \mathbf{M}_L \mathbf{M}_U\]Where \(\mathbf{M}_L\) is lower triangular, and \(\mathbf{M}_U\) is upper triangular. So our equation now reads:
\[\mathbf{LR} = \mathbf{M}_L\mathbf{M}_U\]Triangular matrices are super easy to invert, so invert the right hand side to get:
\[\mathbf{M}_L^{-1}\mathbf{LR}\mathbf{M}_U^{-1} = \mathbf{I}\]Now, since we want left and right eigenvectors that are biorthonormal, we can replace the identity:
\[\mathbf{M}_L^{-1}\mathbf{LR}\mathbf{M}_U^{-1} = \mathbf{L}'\mathbf{R}'\]where the prime indicates our new biorthonormal left and right eigenvectors.
This suggests our new biorthonormal left and right eigenvectors take the form:
\[\mathbf{L}' = \mathbf{M}_L^{-1}\mathbf{L}\]and
\[\mathbf{R}' = \mathbf{R}\mathbf{M}_U^{-1}\]And there you have it! An easy to implement method of biorthonormalizing your eigenvectors. All of the steps have a corresponding LAPACK routine. I’ve found it to be pretty robust. You can go a step further and prove that the new eigenvectors are still eigenvectors of the Hamiltonian. Just plug them in to the eigenvalue equation and you’ll see.
While I think this doesn’t scale any worse than your diagonalization in the first place, it is still a super useful trick. For example, it even works for set of eigenvectors that aren’t full rank (e.g. rectangular). Because you do the LU on the inner product of the left and right eigenvectors, you’ll get a much smaller square matrix the dimension of the number of eigenvectors you have.
Another use you might find with this is if the eigenvectors are nearly biorthonormal (which often happens when you have degenerate eigenvalues). You can do the same trick, but on the subspace of the eigenvectors corresponding to the degenerate eigenvalues. So if you have three degenerate eigenvalues, you can do an LU decomposition plus inversion on a 3x3 matrix.
The end of the year is rapidly approaching, and I’d like to draw your attention to three papers of mine that have been recently published online.
This first paper was a fun collaboration between our theoretical work in the Li group and the experimental work done by Dan Gamelin at the UW. The paper looks at calculated UV-Vis spectra of photodoped and Aluminum-doped ZnO quantum dots. Both types of doping added “extra” electrons to the conduction band, and this resulted in some interesting properties that you can see in both theory and experiment.
What I was most surprised to see was that the “extra” electron simply acted like an electron moving in a spherical potential, just like a hydrogen atom! On the left we have the HOMO, which looks like an s-orbital. Then the two on the right are different LUMOs, one that looks like a p-orbital and another that looks like a d-orbital. We called the “super-orbitals”, since they look like atomic orbitals but exist in these tiny crystal structures. I think these are the first images of DFT-computed “superorbitals” ever published, though the “artificial atom” paradigm has been around for some time. It’s interesting to see complicated systems behaving like the simple quantum mechanical models we study when we first learn quantum mechanics!
J. J. Goings, A. Schimpf, J. W. May, R. Johns. D. R. Gamelin, X. Li, “ Theoretical Characterization of Conduction-Band Electrons in Photodoped and Aluminum-Doped Zinc Oxide (AZO) Quantum Dots,” J. Phys. Chem. C, 2014 , 118, 26584.
This second paper details the implementation of several low-scaling methods for computing excited states (e.g. computing UV/Vis absorption). The idea was to take the highly accurate EOM-CCSD equations and simplify them using perturbation theory. That way, we could keep most of the accuracy of the method, while still having fast enough equations that large molecules could be studied. The resulting equations are closely related to other methods like CC2, CIS(D), and ADC(2), and I showed how they all relate to each other. We compared the performance to EOM-CCSD as well as experimental values. The results were promising, and the perturbative equations performed particularly well for Rydberg states. CC2, on the other hand, performs great for valence excitations.
J. J. Goings, M. Caricato, M. Frisch, X. Li, “ Assessment of Low-scaling Approximations to the Equation of Motion Coupled-Cluster Singles and Doubles Equations,” J. Chem. Phys., 2014 , 141, 164116.
Ab Initio Non-Relativistic Spin Dynamics
Finally, in this paper, we extended the generalized Hartree-Fock method to the time domain. In this proof-of-concept paper, we showed how a magnetic field can guide the spin dynamics of simple spin-frustrated systems. The key is reformulating the real-time time-dependent Hartree-Fock equations in the complex spinor basis. This allows the spin magnetizations on each atom to vary as a response to, say, an externally applied magnetic field. Here’s an example with a lithium trimer. Initially (left picture), all the spin magnetizations (that is, the spatial average of spin magnetic moment) point away from each other. Then, applying a static magnetic field (right picture) into the plane of the molecule causes each magnetization to precess at the Larmor frequency. The precession is shown in picoseconds by the colorization.
It’s really a beautiful idea in my opinion, and there is so much more to be done with it. For example, in our simple ab initio model, the spins only “talk” through Pauli repulsion, so they behave more or less independently. What would happen if we include spin-orbit coupling and other perturbations? That remains to be seen.
F. Ding, J. J. Goings, M. Frisch, X. Li, “ Ab Initio Non-Relativistic Spin Dynamics,” J. Chem. Phys., 2014 , 141, 214111.