# Symmetry analysis for atom arrangements in a crystal

For a while now I've been meaning to post some of my Haskell code here so that someone can tell me what parts of the language/base library I've been completely overlooking. This is the first thing I've brought into a "working/finished" state that isn't a PE problem, and unfortunately it is a tad big.

It is rather mathematical in nature, and contains multiple nontrivial alglorithms which could almost certainly be improved... but I'm hoping to focus less on the algorithms and more on the minutae. I.e. how is my usage of Haskell?

What does it do? Some symmetry analysis for materials research. If you run it, it will print out all of the possible supercell shapes containing 8 16 atoms which are unique under symmetry for a diamond cubic crystal.

Last minute Hawthorne-effect edits aside, the majority of this code was written with the expectation that I would be the only person to ever look at it. I can only hope that at least one soul is brave enough to continue after reading that statement...

Directory structure

$tree -P '*.hs' . ├── Main.hs └── My ├── Common.hs ├── GroupTheory.hs ├── IntegerRref.hs └── Matrix.hs  Main.hs Mostly functions specific to the problem. Here you'll see from all the commented code that my workflow for debugging largely centers around modifying main and recompiling. {-# OPTIONS_GHC -fno-ignore-asserts #-} {-# LANGUAGE BangPatterns #-} import qualified Math.NumberTheory.Primes.Factorisation as Factor -- package arithmoi import qualified Data.Set as Set -- package containers import Data.Set (Set) import qualified Data.List as List import Control.Exception import Text.Printf import Debug.Trace import My.Matrix import My.GroupTheory import My.IntegerRref import My.Common(decorate) ------------------------------------- -- Factorization -- positive, ordered tuples (a,b) such that a*b == x factorPairs :: Integer -> [(Integer, Integer)] factorPairs x = decorate (x div)$ Set.toList $Factor.divisors x -- ordered n-tuplets of positive factors fs such that product fs == x factorTuplets :: Int -> Integer -> [[Integer]] factorTuplets n x = (assert (x>0))$ -- true for expected inputs in this program
case n compare 1 of
LT -> error "factorTuplets: n < 1"
EQ -> [[x]]
GT -> concatMap forPair (factorPairs x)
where forPair (a,b) = map (a:) $factorTuplets (pred n) b ------------------------------------- -- Groups and sets universeForDiagonal :: Vec -> [Mat] universeForDiagonal [] = [[]] universeForDiagonal (a:bcs) = assert (a>0)$ do
let left = a:(fmap (const 0) bcs)
inner <- universeForDiagonal bcs
top <- sequence $fmap (\b -> [0..b-1]) bcs return$ (prependCol left.prependRow top) inner

universeForVolume :: Integer -> [Mat]
universeForVolume vol = concatMap universeForDiagonal (factorTuplets 3 vol)

-------------------------------------
-- Supercell symmetry group

scGenerators :: [Mat]
scGenerators =
-- Twofold rotations
[[0,1,0],[1,0,0],[-1,-1,-1]]:
[[-1,0,0],[0,-1,0],[1,1,1]]:
[[1,0,0],[0,0,1],[-1,-1,-1]]:
-- Threefold rotation
[[0,0,1],[1,0,0],[0,1,0]]:
[]

scMul :: Mat -> Mat -> Mat
scMul = mulMatMat

-- sc group matrices are written to operate on a matrix whose columns are the sc vecs,
-- but we have them in columns;
scAction :: Mat -> Mat -> Mat
scAction g x = integerRref $mulMatMat x (List.transpose g) scGroup :: Set Mat scGroup = generateGroup scMul scGenerators ------------------------------------- -- Random nonsense assertEq :: (Eq a, Show a) => a -> a -> b -> b assertEq expected actual | expected /= actual = error (printf "\n Expected: %s \n Actual: %s" (show expected) (show actual)) | otherwise = id tests :: () tests = id .(assertEq 155$ length $universeForVolume 8) .(assertEq [[-10,5],[2,3]]$ opAddMultiple 3 1 0 $[[-16,-4],[2,3]]) .(assertEq 77$ length $equivalenceClasses scAction (Set.toList scGroup) (universeForVolume 20))$()

main :: IO ()
main = do
let !_ = tests
--  mapM_ print $factorTuplets 3 10 -- mapM_ print$ upperTriangleFromDiags [[1,2,3,4],[5,6,7],[8,9],[1]]
--  mapM_ print $concat$ universeForDiagonal [1,2,3]
--  mapM_ print $opAddMultiple (-2) 2 1$ [[1,4,7,6],[2,3,1,6],[7,7,7,2],[1,2,1,1]]
--  mapM_ print $fmap (getDiag [[1,4,7,6],[2,3,1,6],[7,7,7,2],[1,2,1,1]])$ [-3..3]
--  mapM_ print $fmap (getDiag [[1,4,7,6],[2,3,1,6],[7,7,7,2],[1,2,1,1]])$ [-3..3]
--  mapM_ print $integerRref [[1,1,1],[5,2,2],[3,3,4]] -- mapM_ print$ integerRref [[5,2,2],[1,1,1],[3,3,4]]
--  print $Set.size scGroup -- mapM_ print$ equivalenceClasses scAction (Set.toList scGroup) (universeForVolume 77)
--  mapM_ print $integerRref [[0,1,0],[0,0,1],[4,0,0]] mapM_ print$ map Set.findMin $equivalenceClasses scAction (Set.toList scGroup) (universeForVolume 8)  My/GroupTheory.hs A small number of algorithms related to group theory. module My.GroupTheory ( generateGroup, equivalenceClasses, checkedEquivalenceClasses, ) where import qualified Data.Set as Set -- package containers import Data.Set (Set) setUnionAll :: Ord a => [Set a] -> Set a setUnionAll = foldl Set.union Set.empty -- Generate all elements of a finite group from a generating subset. generateGroup :: Ord g => (g -> g -> g) -> [g] -> Set g generateGroup _ [] = error "generateGroup: empty group" generateGroup mul generators = loop (Set.fromList generators) Set.empty where loop recent output | Set.null recent = output | otherwise = loop new (output Set.union new) where new = (Set.difference output)$ Set.fromList $fmap (\[a,b] -> a mul b)$
sequence [Set.toList recent, generators]

equivalenceClasses :: Ord x => (g -> x -> x) -> [g] -> [x] -> [Set x]
equivalenceClasses _ [] _  = error "equivalenceClasses: empty group"
equivalenceClasses action group xs =
loop (Set.fromList xs)
where loop remaining
| Set.null remaining = []
| otherwise = newClass:loop (remaining Set.difference newClass)
where newClass = Set.fromList $fmap (action (Set.findMin remaining)) group checkedEquivalenceClasses :: Ord x => (g -> x -> x) -> [g] -> [x] -> [Set x] checkedEquivalenceClasses action group xs = validate classes where classes = equivalenceClasses action group xs validate = validateClosed.validateDisjoint -- each element of xs should appear once and only once in the classes validateDisjoint = case totalOutCount compare uniqueOutCount of LT -> error "checkedEquivalenceClasses: internal error" GT -> error "checkedEquivalenceClasses: classes not disjoint" EQ -> id -- the union of the classes must equal xs validateClosed | (not . Set.null) (uniqueIn Set.difference uniqueOut) = error "checkedEquivalenceClasses: internal error in equivalenceClasses" | (not . Set.null) (uniqueOut Set.difference uniqueIn) = error "checkedEquivalenceClasses: group action not closed on xs" | otherwise = id uniqueIn = Set.fromList xs uniqueOut = setUnionAll classes uniqueOutCount = Set.size uniqueOut totalOutCount = sum (map Set.size classes) -- findInverses :: Set g -> Map g g -- findIdentity :: Set g -> g -- validateClosure :: Set g -> Set g -- O(n^2) -- validateAssociativity :: Set g -> Set g -- O(n^3) -- validateGroup :: Set g -> Set g -- validateAction :: Set g -> [x] -> ()  My/Matrix.hs Defines operations on vectors and matrices, "implemented" as little more than lists. module My.Matrix ( Mat, Vec, prependCol, prependRow, deleteCol, overRow, overLowerRight, getCol, getDiag, upperTriangleFromDiags, height, width, innerProd, mulMatMat, mulMatVec, ) where import qualified Data.List as List import Control.Exception import My.Common(listSet,deleteAt,zipWithExact) type Mat = [[Integer]] type Vec = [Integer] prependCol :: Vec -> Mat -> Mat prependRow :: Vec -> Mat -> Mat prependCol col mat = zipWith (:) col mat prependRow = (:) deleteCol :: Int -> Mat -> Mat deleteCol i = fmap (deleteAt i) -- think "over (_!! i)", if _!! were some sort of Lens for lists overRow :: Int -> (Vec -> Vec) -> Mat -> Mat overRow i f mat = listSet i (f (mat!!i)) mat -- applies a function to the lower right submatrix excluding nDropped -- rows and columns. overLowerRight :: Int -> (Mat -> Mat) -> Mat -> Mat overLowerRight nDropped f mat = let (top, bottom) = List.splitAt nDropped mat in let (botLs, botRs) = unzip$ map (List.splitAt nDropped) bottom in
top ++ (zipWith (++) botLs (f botRs))

getCol :: Int -> Mat -> Vec
getCol i = fmap (!!i)

getDiag :: Int -> Mat -> Vec
getDiag n rows = case n compare 0 of
LT -> getDiag 0 (drop (-n) rows)
GT -> getDiag 0 (fmap (drop n) rows)
EQ -> fmap (\(i,row) -> row!!i) $take w (zip [0,1..] rows) where w = (length.head) rows upperTriangleFromDiags :: [Vec] -> Mat upperTriangleFromDiags [] = [] upperTriangleFromDiags diags = assert checks$ topRow:otherRows where
checks = length diags == length (head diags)
topRow = map head diags
otherRows = map (0:) $upperTriangleFromDiags (map tail (init diags)) height :: Mat -> Int height = length width :: Mat -> Int -- We do not explicitly store a width, so none can be determined from -- a zero row matrix. That said, I never plan to use one, so better -- safe than sorry: width [] = error "width: null matrix" width mat = (length.head) mat innerProd :: Vec -> Vec -> Integer innerProd a b = sum$ zipWithExact (*) a b

mulMatMat :: Mat -> Mat -> Mat
mulMatMat a b
| (width a) /= (height b) = error "mulMatMat: dimension"
| otherwise = [[innerProd row col | col<-List.transpose b] | row <- a]

mulMatVec :: Mat -> Vec -> Vec
mulMatVec m v
| (width m) /= (length v) = error "mulMatVec: dimension"
| otherwise = map (innerProd v) m


My/IntegerRref.hs

The largest piece of code, implementing an algorithm based on various handwritten proofs. It is very closely related to the Hermite normal form of a matrix.

module My.IntegerRref (
opAddRow, opSubRow, opNegate2, opAddMultiple,
integerRref,
validateIrref,
) where

import qualified Data.List as List
import Control.Exception

import My.Matrix
import My.Common(listSet,pDiv,compose)

---------------------------------------------------
-- This module implements an analogue to Reduced Row Echelon Form where
-- the only primitive symmetry operations are the following:
--   * adding one row into another, different row.
--   * subtracting one row from another, different row.
--
-- These operations preserve the determinant of a matrix.

opSubRow :: Int -> Int -> Mat -> Mat
opSubRow = opAddMultiple (-1)
opAddRow :: Int -> Int -> Mat -> Mat
opAddRow = opAddMultiple 1

-- We can't negate one row, but we can negate two:
opNegate2 :: Int -> Int -> Mat -> Mat
opNegate2 i1 i2 = overRow i1 (fmap negate).overRow i2 (fmap negate)

opAddMultiple :: Integer -> Int -> Int -> Mat -> Mat
opAddMultiple b src dest mat
| dest == src = error "rowAdd: src == dest (not volume-conserving)"
| otherwise = listSet dest newRow mat where
newRow = zipWith (+) (mat!!dest) (fmap (b*) (mat!!src))

---------------------------------------------------

-- For a square, invertible, integer matrix, produces the unique matrix that is:
--  * Reachable from the original by a finite sequence of operations consisting
--    of adding an integer multiple of one row into another (different) row.
--  * Is upper triangular.
--  * Is entirely nonnegative with the SINGLE possible exception of the
--    lower-right most element.
--  * For each column, the absolute value of the element on the main diagonal
--    is strictly greater than all other values in the column.
integerRref :: Mat -> Mat
integerRref [] = []
integerRref mat = validateIrref $compose operations$ partialRref mat
where operations = do
rPivot <- [0..(length mat)-1]
rMod <- [0..rPivot-1]
return $(reduceRowModuloRow rPivot rPivot) rMod -- This produces the correct values in the lower triangular part of the matrix -- I.e. the main diagonal is positive with the possible exception of the last -- element, and off-diagonals below the main diagonal are zero. partialRref :: Mat -> Mat partialRref = rec where rec [] = [] rec mat = (fixRest.fixLeft) mat where -- change column into [g,g,g..], and then into [g,0,0..] fixLeft = compose$
(reduceColumnToConstant 0):
map (opSubRow 0) [1..(length mat)-1]
-- recurse
fixRest = overLowerRight 1 rec

-- reduces a column [a,b,c..] into [g,g,g..], where g = gcd (a,b,c..).
--  by performing primitive row ops.  As long as the original matrix
--  is larger than 1x1, g will be non-negative.
reduceColumnToConstant :: Int -> Mat -> Mat
reduceColumnToConstant _ [] = error "reduceColumnToConstant: empty mat" -- no columns!
reduceColumnToConstant _ [row] = [row] -- do NOT call ensureColumnNonNegative
reduceColumnToConstant c mat'@(_:_) = (verify.compose operations) mat' where
verify mat = assert (all (==mat!!0!!c) (getCol c mat)) mat

-- The column is made positive to simplify working with the gcd.
-- Then two passes are made where each row is reduced with the top row by
--   adding multiples of one to another until they share a value in column c.
-- The first pass puts the final correct value in the first and last rows;
--   the second pass propagates it to the other rows.
operations = (ensureColumnNonNegative c):(onePass ++ onePass) where
onePass = map (reduceTwoRows 0) [1..(length mat')-1]

reduceTwoRows r1 r2 = verify.loop where
-- like a bizarre variant of euclid's algorithm, where we're performing
--  operations on entire rows, and we're not allowed to swap their order
verify mat = assert (mat!!r1!!c == mat!!r2!!c) mat
loop mat = case (mat!!r1!!c, mat!!r2!!c) of
(0, 0) -> mat
(0, _) -> (opAddRow r2) r1 mat
(_, 0) -> (opAddRow r1) r2 mat
(g, h) -> case (g compare h) of
-- use lesser to reduce larger
LT -> loop $(reduceRowModuloRow c r1) r2 mat GT -> loop$ (reduceRowModuloRow c r2) r1 mat
EQ -> mat

ensureColumnNonNegative :: Int -> Mat -> Mat
ensureColumnNonNegative c = verify.loop where
verify mat = assert (all (0<=) (getCol c mat)) mat
loop mat = case (positives, negatives, zeros) of

-- Column is non-negative
(_,  [],   _) -> mat

-- Column is hopeless
([], [_], []) -> error "ensureColumnPositive: only one row and value is negative"

-- At least one negative and one other: we can negate them
([], n1:n2:_, _) -> loop $opNegate2 n1 n2 mat ([], [n1], z1:_) -> loop$ opNegate2 n1 z1 mat

-- At least one positive p:
--    use it to bring the negatives into the range 0..p-1
(p:_, ns@(_:_), _) -> compose (fmap fixRow ns) mat
where fixRow = reduceRowModuloRow c p

where
negatives = List.findIndices ( <0) column
zeros     = List.findIndices (==0) column
positives = List.findIndices ( >0) column
column = getCol c mat

-- add a multiple of row pivotRow into row modRow that reduces
--  the element in column col into the range [0..|pivot| - 1].
-- The pivot must be nonzero.
reduceRowModuloRow :: Int -> Int -> Int -> Mat -> Mat
reduceRowModuloRow col pivotRow modRow mat =
let moderand = mat!!modRow!!col in
let pivot = mat!!pivotRow!!col in
let d = moderand pDiv pivot in
opAddMultiple (-d) pivotRow modRow mat

-- note: this assumes the matrix is square and invertible, so that
--       the pivots MUST be on the main diagonal
-- warning: this inspects every element, causing any deferred lazy
--          computations to occur conditionally based on whether
--          or not assertions are enabled.  I imagine there is the
--          potential for this to hide nasty memory bugs...
validateIrref :: Mat -> Mat
validateIrref mat = checkDims $checkDiag$ checkZeros $checkReduced$ mat
where
checkDims = assert $all ((==length mat).length) mat -- the final pivot is the only element allowed to be negative in the -- entire matrix checkDiag = assert$ (all (0<) (init pivots)) && (0 /= last pivots)

checkZeros = assert $all test$ zip [0,1..] mat
where test (i,row) =
i == (length $takeWhile (0==)$ row)

checkReduced = assert $all test$ List.zip3 [0,1..] columns pivots
where test (i,col,pivot) =
all (\x -> 0 <= x && x < abs pivot) $(take i col) columns = List.transpose mat pivots = getDiag 0 mat  My/Common.hs Various helper functions that I have deemed "missing" from the Haskell base library. (Your job is to tell me why they aren't!) module My.Common ( deleteAt, listSet, zipWithExact, decorate, compose, indices, windows, traceWith, pDivMod,pDiv,pMod, ) where import Debug.Trace import qualified Data.List as List deleteAt :: Int -> [a] -> [a] deleteAt _ [] = error "deleteAt: i >= length" deleteAt i xs@(_:xt) | i < 0 = error "deleteAt: i < 0" | i == 0 = xt | otherwise = deleteAt (i-1) xs listSet :: Int -> a -> [a] -> [a] listSet _ _ [] = error "listSet: index (or empty list)" listSet 0 new (_:xs) = new:xs listSet n new (x:xs) = x:listSet (n-1) new xs zipWithExact :: (a -> b -> c) -> [a] -> [b] -> [c] zipWithExact f as' bs' = iter as' bs' where iter [] [] = [] iter [] (_:_) = error "zipWithExact: first list ended early" iter (_:_) [] = error "zipWithExact: second list ended early" iter (a:as) (b:bs) = (f a b):iter as bs decorate :: Functor f => (a -> b) -> f a -> f (a,b) decorate f = fmap (\x -> (x, f x)) -- compose [f,g,h..] x == (..h.g.f) x compose :: [(a -> a)] -> a -> a compose = foldl (flip (.)) id -- equivalent to [0..(length xs-1)] but possibly less painful? indices :: [a] -> [Int] indices xs = map fst$ zip [0..] xs

-- overlapping windows
windows :: Int -> [a] -> [[a]]
windows n xs = takeWhile ((n==).length) $fmap (take n)$ List.tails xs

-- threads a value through a function, printing the output of the
--  function and returning the original value
traceWith :: (a -> String) -> a -> a
traceWith f x = trace (f x) x

-- HACK; Debug.Trace apparently should have this but mine doesn't.
-- (my base library is probably still the one from the Canonical repos >_>)
traceShowId :: Show a => a -> a
traceShowId = traceWith show

-- another variant of divMod and quotRem satisfying the property
-- that pMod is never negative, even if the divisor is.
pDivMod :: Integral a => a -> a -> (a,a)
pDivMod a b = (d, a-b*d)
where d = (signum b) * (a div (abs b))

pDiv :: Integral a => a -> a -> a
pMod :: Integral a => a -> a -> a
pDiv a b = fst $pDivMod a b pMod a b = snd$ pDivMod a b


Specific points of concern:

• In IntegerRef you'll see the use of a helper function called compose to string long chains of operations together into one function that performs them in sequence. I understand that this is precisely the sort of problem that Monads are meant to solve, but I'm not sure how one could help. (the only ones I really understand are [] and Maybe)
• Something I mention in the notes above validateIrref (in IntegerRref.hs); a typical debugging strategy of mine is to verify post-conditions of functions through expensive checks which are only enabled in debug mode. But in Haskell, with lazy evaluation, that might cause the "strictness" of a function to depend on whether assertions are enabled! Seems troubling...
• Better ways to write more things in a pointfree style. Or conversely, places in my code where I used a pointfree style to the detriment of readability.
• Any "anti-idioms" I use that I should be aware of
• Why not use the aptly named matrix library? 1) it's much faster, 2) you shouldn't need to worry whether their functions are correct so can debug less. Jun 9 '16 at 15:05
• @MichaelKlein a valid question. Obvious true answer aside (which is "I didn't look for one"), I think part of it is that most of my proofs were structured recursively and so lists seemed a natural fit. The matrix-related type aliases and methods are largely things I pulled out as an afterthought when it began to bother me that code working on columns tended to look so different from code working on rows. Jun 10 '16 at 0:15

## 1 Answer

Here's a little to start with:

Safe.Exact implements zipWithExact.

Lens implements listSet n as set (ix n), and overRow n as over (ix n). These don't error on being out of range, instead doing nothing. listSet is a bad name because list is also a verb and set is also a noun.

overLowerRight n is over (foldr (.) id (replicate n _tail) . each . foldr (.) id (replicate n _tail) . each).

generateGroup should be called generateAbelianMagma or generateSemigroup, because you aren't generating the neutral element or inverses, and are using either commutativity or associativity to only append the original generators to any new elements, and only to the right.

generateSemigroup :: Ord g => (g -> g -> g) -> [g] -> Set g
generateSemigroup op generators = foldr foo Set.empty generators where
foo :: g -> Set g -> Set g
foo x set = if S.member x set
then set
else foldr foo (insert x set) $(op x) <$> generators

• Hmm... Perhaps I should have called it generateFiniteGroup? I think that ought to at least address the inverses and identity problem because then each generator has the identity in its cyclic subgroup. I'm still trying to digest your new definition, but it's impressive! Jun 6 '16 at 17:54
• Ah, that explains why you error out on an empty generator list - the neutral element isn't generated then. And for infinite groups, this doesn't halt anyway, I guess. Jun 7 '16 at 21:30