The code implements a specific equation in the process of calculating the elastic constants of an arbitrary 3 dimensional cell from the fluctuations in the cell. The cell definition is contained in a 3x3 matrix. I read in a data file containing many lines, where each line holds the necessary information to construct the cell matrix (The cell matrix is upper triangular, so it’s only 6 data points). The average cell matrix is calculated (as you’d expect, more or less), then there’s a matrix equation which defines some necessary quantities for the calculation of the relative fluctuation of the matrices. At that point I calculate the strain on the matrix which is an outer matrix product (and hence promotes me from using a 3x3 matrix - a rank 2 tensor to a 3x3x3x3 4th rank tensor). The code I’ve posted is simply the contraction of that 4th rank tensor (due to symmetry) into a more amenable form. It turns out that due to the symmetry of real world physical systems, the 4th rank tensor can be contracted into a 6x6 matrix, which is naturally composed, due to symmetries, of 4 3x3 matrices (I excised the calculations of cInd1, cInd2, rInd1, and rInd2, but they have symmetric forms which index into the 4th rank tensor as part of this contraction). You can see some of this symmetry in the layout of the variables, especially if you know that the form is:
X(6x6) = |M1(3x3) M2(3x3)|
|M3(3x3) M4(3x3)|
So the loops I’ve posted are basically constructing the M1, M2, M3 and M4 matrices via the appropriate tensor contraction (which consists of index permutations to generate cInd1, cInd2, rInd1, and rInd2, and in some case addition operators).
eToS is the conversion factor from strain (e) to compliance (S) – i.e. the dimensional constants. (Of course in trying to be more succinct I cut out the calculation of eToS where I had variables defined like kdBoltzmannSGI which is a double precision constant that defines Boltzmann’s constant in SGI units, using such obtusely named variables as dVolAvg for the average volume of the system. ;) )
Much of it likely stems from my primarily using Python now; I’m not even sure I understand some of the syntax you’re using anymore! I don’t think I’ve even written an old fashioned C-stye for-loop in years…
Perhaps. I suspect much of it is that I’ve named variables for a very specific set of knowledge which you’d neither guess nor likely be familiar with if you did guess. It comes from converting standard notation equations directly to code. I could use longer descriptive names, but frankly I shudder at the thought of N column wide code like so:
ReducedComplianceMatrix_UpperLeftSubmatrix[CurrentWindowIndex](Row,Column)=ComplianceTensor(TensorIndexI,TensorIndexJ,ReducedRowIndex1,ReducedRowIndex2)+ComplianceTensor(TensorIndexI,TensorIndexJ,ReducedRowIndex2,ReducedRowIndex1);
Is something like that really easier to read? Is there a better intermediate that optimizes readability? It has more characters, but is it magically more meaningful if you don’t know what a ReducedComplianceMatrix is or what a ReducedRowIndex1 is? (In fact the latter seems more confusing to me; how do you define the “row” in a 4 dimensional data structure?)
I’m not trying to be argumentative or throw out tons of minutiae here to push you off; I’m genuinely curious if this information helps, because if it does I’ll happily code in a style that makes things more understandable. I just never really considered the thought that this stuff would be understandable anyway if you weren’t sitting with a paper reference in one hand to compare to, and at that point my variables are named pretty concurrently with the literature source they’re taken from.
