module Crystallography.HallSymbols (
fromHallSymbols,
fromHallSymbols',
generatorsOfHallSymbols,
hallSymbols,
parser,
generators,
LatticeSymbol,
MatrixSymbol,
OriginShift,
) where
import Data.Maybe
import Data.Char (isSpace)
import Data.List (nub,sort,elemIndex)
import Data.Fixed (mod')
import Data.Ratio
import Data.Matrix hiding (identity,matrix,(<|>))
import qualified Data.Matrix as M (identity)
import Text.ParserCombinators.Parsec
import Text.ParserCombinators.Parsec.Error
type LatticeSymbol = (Bool,Char)
type NFold = Int
data MatrixSymbol
= MatrixSymbol {
minusSign :: Bool,
nFoldBody :: NFold,
nFoldSubscript :: Maybe Int,
nFoldDiagonal :: Maybe Char,
axisOfRotation :: Maybe Char,
translationVector :: String
} deriving Show
type OriginShift = (Integer,Integer,Integer)
latticeSymbol :: CharParser () LatticeSymbol
latticeSymbol = do
sign <- optionMaybe (oneOf "-")
b <- oneOf "PABCIRF"
return (isJust sign,b)
matrixSymbol :: CharParser () MatrixSymbol
matrixSymbol = do
sign <- optionMaybe (oneOf "-")
b <- oneOf "12346"
s <- optionMaybe (oneOf "12345")
d <- optionMaybe (oneOf "'*\"")
a <- optionMaybe (oneOf "xyz")
t <- vectorSymbol
if isJust s && isJust d
then
fail "Symbols that can not coexist."
else
return $ MatrixSymbol (isJust sign) (read [b]) (read . (: []) <$> s) d a t
vectorSymbol :: CharParser () String
vectorSymbol = many (oneOf "abcnuvwd")
originShift :: CharParser () OriginShift
originShift = do
string "("
va <- integer
space
vb <- integer
space
vc <- integer
string ")"
return (va,vb,vc)
integer :: CharParser () Integer
integer = signed <|> unsigned
where
digits = many1 digit
signed = do
a <- char '-'
num <- digits
return (read $ a : num)
unsigned = do
num <- digits
return (read num)
space' :: CharParser () Char
space' = oneOf " _"
parser :: CharParser () ( LatticeSymbol, [MatrixSymbol], OriginShift )
parser = do
l <- latticeSymbol
nat <- many1 (try ms)
v <- option (0,0,0) os
return (l, nat, v)
where
ms = do
space'
matrixSymbol
os = do
space'
originShift
generators :: CharParser () [Matrix Rational]
generators = do
raw <- parser
return $ decodeSymbols raw
hallSymbols :: CharParser () [Matrix Rational]
hallSymbols = do
g <- generators
let equivalentPositions = generate g
if null equivalentPositions
then
fail "something happen when decode or generate process."
else
return equivalentPositions
generateLikeITC :: CharParser () [Matrix Rational]
generateLikeITC = do
raw <- parser
let equivalentPositions = likeITC raw
if null equivalentPositions
then
fail "something happen when decode or generate process."
else
return equivalentPositions
likeITC :: ( LatticeSymbol, [MatrixSymbol], OriginShift ) -> [Matrix Rational]
likeITC raw = if null h then [] else r
where
constructMatrices' (l,nat,v) = fmap (mapOriginShfit v) [lattice l, map matrix nat]
[a,b] = constructMatrices' . restoreDefaultAxis $ raw
h = generate $ [identity] ++ b
r = generate $ h ++ a
fromHallSymbols :: String -> Either ParseError [Matrix Rational]
fromHallSymbols s = parse generateLikeITC ("while reading " ++ show s) s
fromHallSymbols' :: String -> [Matrix Rational]
fromHallSymbols' s = case fromHallSymbols s of
Left e -> error $ show e
Right mm -> mm
generatorsOfHallSymbols :: String -> [Matrix Rational]
generatorsOfHallSymbols s = case gg s of
Left e -> error $ show e
Right mm -> mm
where
gg s = parse generators ("while reading " ++ show s) s
fromHallSymbols''' :: String -> Either ParseError [Matrix Rational]
fromHallSymbols''' s = parse hallSymbols ("while reading " ++ show s) s
decodeSymbols :: ( LatticeSymbol, [MatrixSymbol], OriginShift ) -> [Matrix Rational]
decodeSymbols = constructMatrices . restoreDefaultAxis
constructMatrices :: (LatticeSymbol, [MatrixSymbol], OriginShift) -> [Matrix Rational]
constructMatrices (l,nat,v) = mapOriginShfit v $ lattice l ++ map matrix nat
restoreDefaultAxis :: (LatticeSymbol, [MatrixSymbol], OriginShift) -> (LatticeSymbol, [MatrixSymbol], OriginShift)
restoreDefaultAxis ( l, nat, v ) = ( l, mapDA nat, v )
mapDA :: [MatrixSymbol] -> [MatrixSymbol]
mapDA = mapDA' 1 0
mapDA' :: Int -> NFold -> [MatrixSymbol] -> [MatrixSymbol]
mapDA' _ _ [] = []
mapDA' o p (x:xs) = da o p x : mapDA' (succ o) (nFoldBody x) xs
da :: Int -> NFold -> MatrixSymbol -> MatrixSymbol
da 0 _ _ = error "order parameter must be >1."
da 1 preceded (MatrixSymbol s n1 n2 n3 Nothing v)
= MatrixSymbol s n1 n2 n3 (Just 'z') v
da 2 preceded (MatrixSymbol s 2 Nothing Nothing Nothing v)
| preceded `elem` [2,4] = MatrixSymbol s 2 Nothing Nothing (Just 'x') v
| preceded `elem` [3,6] = MatrixSymbol s 2 Nothing (Just '\'') (Just 'z') v
da 3 _ (MatrixSymbol s 3 Nothing Nothing a v)
= MatrixSymbol s 3 Nothing (Just '*') a v
da order preceded (MatrixSymbol s n1 n2 n3 Nothing v)
= MatrixSymbol s n1 n2 n3 (Just 'z') v
da order preceded (matrixSymbol@(MatrixSymbol _ n1 _ _ _ _))
= matrixSymbol
generate :: [Matrix Rational] -> [Matrix Rational]
generate mm = gn 0 mm mm
gn :: Int -> [Matrix Rational] -> [Matrix Rational] -> [Matrix Rational]
gn n s m | length (nub m) == length mm = mm
| n > 10 = []
| otherwise = gn (succ n) (nub s) mm
where
mm = nub . map (modulus1 . foldl1 multStd) . sequenceA $ [s,m]
setTransform :: Matrix Rational -> [Rational] -> Matrix Rational
setTransform m l = let t = fromList 3 1 l
(a,b,c,d) = splitBlocks 3 3 m
in joinBlocks (a,t,c,d)
mapOriginShfit :: (Integer,Integer,Integer) -> [Matrix Rational] -> [Matrix Rational]
mapOriginShfit (va,vb,vc) = map (sn va vb vc)
where
sn 0 0 0 m = m
sn va vb vc m = let v = [ va % 12, vb % 12, vc % 12 ]
n = setTransform identity v
o = setTransform identity (fmap negate v)
in modulus1 $ foldl1 multStd [n,m,o]
modulus1 :: Real a => Matrix a -> Matrix a
modulus1 m = let (a,b,c,d) = splitBlocks 3 3 m
in joinBlocks (a,fmap (`mod'` 1) b,c,d)
ng33 :: Matrix Rational -> Matrix Rational
ng33 m = let (a,b,c,d) = splitBlocks 3 3 m
in joinBlocks (negate a,b,c,d)
identity :: Matrix Rational
identity = M.identity 4
mapNg33 :: [Matrix Rational] -> [Matrix Rational]
mapNg33 l = nub [ g m | g <- [ id, ng33 ], m <- l ]
centroSymmetric :: [Matrix Rational] -> [Matrix Rational]
centroSymmetric = mapNg33
set :: [[Rational]] -> [Matrix Rational]
set = map (setTransform identity)
lattice :: LatticeSymbol -> [Matrix Rational]
lattice (False, l) = tbl1 l
lattice ( True, l) = centroSymmetric (tbl1 l)
tbl1 :: Char -> [Matrix Rational]
tbl1 'P' = set [[ 0, 0, 0 ]]
tbl1 'A' = set [[ 0, 0, 0 ], [ 0, 1%2, 1%2 ]]
tbl1 'B' = set [[ 0, 0, 0 ], [ 1%2, 0, 1%2 ]]
tbl1 'C' = set [[ 0, 0, 0 ], [ 1%2, 1%2, 0 ]]
tbl1 'I' = set [[ 0, 0, 0 ], [ 1%2, 1%2, 1%2 ]]
tbl1 'R' = set [[ 0, 0, 0 ], [ 1%3, 2%3, 2%3 ], [ 2%3, 1%3, 1%3 ]]
tbl1 'F' = set [[ 0, 0, 0 ], [ 0, 1%2, 1%2 ], [ 1%2, 0, 1%2 ], [ 1%2, 1%2, 0 ]]
tbl1 c = error $ show c
tv :: String -> [Rational]
tv ch | null ch = [0,0,0]
| otherwise = fmap (`mod'` 1) $ foldl1 (zipWith (+)) $ map tbl2 ch
tbl2 :: Char -> [Rational]
tbl2 'a' = [ 1%2, 0, 0 ]
tbl2 'b' = [ 0, 1%2, 0 ]
tbl2 'c' = [ 0, 0, 1%2 ]
tbl2 'n' = [ 1%2, 1%2, 1%2 ]
tbl2 'u' = [ 1%4, 0, 0 ]
tbl2 'v' = [ 0, 1%4, 0 ]
tbl2 'w' = [ 0, 0, 1%4 ]
tbl2 'd' = [ 1%4, 1%4, 1%4 ]
tbl2 c = error $ show c
matrix :: MatrixSymbol -> Matrix Rational
matrix (MatrixSymbol False 3 (Just 1) Nothing (Just 'z') _) = setTransform matrix3z [0,0,1%3]
matrix (MatrixSymbol False 3 (Just 2) Nothing (Just 'z') _) = setTransform matrix3z [0,0,2%3]
matrix (MatrixSymbol False 4 (Just 1) Nothing (Just 'z') _) = setTransform matrix4z [0,0,1%4]
matrix (MatrixSymbol False 4 (Just 3) Nothing (Just 'z') _) = setTransform matrix4z [0,0,3%4]
matrix (MatrixSymbol False 6 (Just 1) Nothing (Just 'z') _) = setTransform matrix6z [0,0,1%6]
matrix (MatrixSymbol False 6 (Just 2) Nothing (Just 'z') _) = setTransform matrix6z [0,0,2%6]
matrix (MatrixSymbol False 6 (Just 4) Nothing (Just 'z') _) = setTransform matrix6z [0,0,4%6]
matrix (MatrixSymbol False 6 (Just 5) Nothing (Just 'z') _) = setTransform matrix6z [0,0,5%6]
matrix (MatrixSymbol s n1 n2 n3 a t) = setTransform (ng33if s $ tbl345 n1 n3 a) (tv t)
ng33if flag = if flag then ng33 else id
matrix3z = tbl345 3 Nothing (Just 'z')
matrix4z = tbl345 4 Nothing (Just 'z')
matrix6z = tbl345 6 Nothing (Just 'z')
tbl345 :: Int -> Maybe Char -> Maybe Char -> Matrix Rational
tbl345 1 Nothing _ = fromLists [[ 1, 0, 0, 0], [ 0, 1, 0, 0], [ 0, 0, 1, 0], [ 0, 0, 0, 1]]
tbl345 2 Nothing (Just 'x') = fromLists [[ 1, 0, 0, 0], [ 0,-1, 0, 0], [ 0, 0,-1, 0], [ 0, 0, 0, 1]]
tbl345 3 Nothing (Just 'x') = fromLists [[ 1, 0, 0, 0], [ 0, 0,-1, 0], [ 0, 1,-1, 0], [ 0, 0, 0, 1]]
tbl345 4 Nothing (Just 'x') = fromLists [[ 1, 0, 0, 0], [ 0, 0,-1, 0], [ 0, 1, 0, 0], [ 0, 0, 0, 1]]
tbl345 6 Nothing (Just 'x') = fromLists [[ 1, 0, 0, 0], [ 0, 1,-1, 0], [ 0, 1, 0, 0], [ 0, 0, 0, 1]]
tbl345 2 Nothing (Just 'y') = fromLists [[-1, 0, 0, 0], [ 0, 1, 0, 0], [ 0, 0,-1, 0], [ 0, 0, 0, 1]]
tbl345 3 Nothing (Just 'y') = fromLists [[-1, 0, 1, 0], [ 0, 1, 0, 0], [-1, 0, 0, 0], [ 0, 0, 0, 1]]
tbl345 4 Nothing (Just 'y') = fromLists [[ 0, 0, 1, 0], [ 0, 1, 0, 0], [-1, 0, 0, 0], [ 0, 0, 0, 1]]
tbl345 6 Nothing (Just 'y') = fromLists [[ 0, 0, 1, 0], [ 0, 1, 0, 0], [-1, 0, 1, 0], [ 0, 0, 0, 1]]
tbl345 2 Nothing (Just 'z') = fromLists [[-1, 0, 0, 0], [ 0,-1, 0, 0], [ 0, 0, 1, 0], [ 0, 0, 0, 1]]
tbl345 3 Nothing (Just 'z') = fromLists [[ 0,-1, 0, 0], [ 1,-1, 0, 0], [ 0, 0, 1, 0], [ 0, 0, 0, 1]]
tbl345 4 Nothing (Just 'z') = fromLists [[ 0,-1, 0, 0], [ 1, 0, 0, 0], [ 0, 0, 1, 0], [ 0, 0, 0, 1]]
tbl345 6 Nothing (Just 'z') = fromLists [[ 1,-1, 0, 0], [ 1, 0, 0, 0], [ 0, 0, 1, 0], [ 0, 0, 0, 1]]
tbl345 2 (Just '\'') (Just 'x') = fromLists [[-1, 0, 0, 0], [ 0, 0,-1, 0], [ 0,-1, 0, 0], [ 0, 0, 0, 1]]
tbl345 2 (Just '"') (Just 'x') = fromLists [[-1, 0, 0, 0], [ 0, 0, 1, 0], [ 0, 1, 0, 0], [ 0, 0, 0, 1]]
tbl345 2 (Just '\'') (Just 'y') = fromLists [[ 0, 0,-1, 0], [ 0,-1, 0, 0], [-1, 0, 0, 0], [ 0, 0, 0, 1]]
tbl345 2 (Just '"') (Just 'y') = fromLists [[ 0, 0, 1, 0], [ 0,-1, 0, 0], [ 1, 0, 0, 0], [ 0, 0, 0, 1]]
tbl345 2 (Just '\'') (Just 'z') = fromLists [[ 0,-1, 0, 0], [-1, 0, 0, 0], [ 0, 0,-1, 0], [ 0, 0, 0, 1]]
tbl345 2 (Just '"') (Just 'z') = fromLists [[ 0, 1, 0, 0], [ 1, 0, 0, 0], [ 0, 0,-1, 0], [ 0, 0, 0, 1]]
tbl345 3 (Just '*') _ = fromLists [[ 0, 0, 1, 0], [ 1, 0, 0, 0], [ 0, 1, 0, 0], [ 0, 0, 0, 1]]
tbl345 a b c = error $ show (a,b,c)