ZEBEDDE (empty) → 0.1.0.0

raw patch · 8 files changed

+314/−0 lines, 8 filesdep +basedep +vectsetup-changed

Dependencies added: base, vect

Files

LICENSE +1/−0
Setup.hs +2/−0
ZEBEDDE.cabal +24/−0
ZEBEDDE/Core/Molecule.lhs +63/−0
ZEBEDDE/Core/ReadXYZ.hs +61/−0
ZEBEDDE/Core/Vector.lhs +81/−0
ZEBEDDE/Quaternion.hs +26/−0
ZEBEDDE/Translation.lhs +56/−0

+ LICENSE view

@@ -0,0 +1,1 @@+all rights reserved Guy Storey & Dewi Lewis

+ Setup.hs view

@@ -0,0 +1,2 @@+import Distribution.Simple
+main = defaultMain

+ ZEBEDDE.cabal view

@@ -0,0 +1,24 @@+-- Initial ZEBEDDE.cabal generated by cabal init.  For further 
+-- documentation, see http://haskell.org/cabal/users-guide/
+
+name:                ZEBEDDE
+version:             0.1.0.0
+synopsis:            Polymer growth simulation method
+-- description:         
+license:             BSD3
+license-file:        LICENSE
+author:              Guy Storey
+maintainer:          guystorey88@hotmail.com
+-- copyright:           
+category:            Math
+build-type:          Simple
+-- extra-source-files:  
+cabal-version:       >=1.10
+
+library
+  exposed-modules:     ZEBEDDE.Translation, ZEBEDDE.Quaternion, ZEBEDDE.Core.Vector, ZEBEDDE.Core.ReadXYZ, ZEBEDDE.Core.Molecule
+  -- other-modules:       
+  other-extensions:    GeneralizedNewtypeDeriving, TypeSynonymInstances, FlexibleInstances
+  build-depends:       base >=4.7 && <4.8, vect >= 0.4.7
+  -- hs-source-dirs:      
+  default-language:    Haskell2010

+ ZEBEDDE/Core/Molecule.lhs view

@@ -0,0 +1,63 @@+
+\section{Atoms and Molecules}
+
+> module ZEBEDDE.Core.Molecule (Atom,
+>                              	element,
+>                               mass,
+>                              	position,
+>                               doToAtom,
+>                               Molecule,
+>                               name,
+>                               basis,
+>                               atoms,
+>                               newMolecule,
+>                               doToMolecule,
+>                               doToMolecule') where
+
+
+> import ZEBEDDE.Core.Vector
+
+\subsection{Atoms}
+
+An atom has a name, a molecular weight and a position. 
+
+> type Atom = (String,Int,Vec)       
+
+Access elements of an atom.
+
+> element    :: Atom -> String
+> element    (a,_,_) = a
+
+> mass     :: Atom -> Int
+> mass     (_,b,_) = b
+
+> position :: Atom -> Vec
+> position (_,_,c) = c
+
+Apply function to the vector component of an atom.
+
+> doToAtom f (string,mass,vec) = (string,mass,f vec)
+
+\subsection{Molecule}
+Create molecule strucure, \\ 3 vecorts specifying location and orientation \\ and it is a name and a list of atoms.
+
+> type Basis = (Vec,Vec,Vec)
+> type Molecule = (String,Basis,[Atom])
+> name  (n,_,_) = n
+> basis (_,b,_) = b
+> atoms (_,_,a) = a
+
+> getBasis :: [Atom] -> Basis
+> getBasis atoms = (position (atoms!!0),position (atoms!!1),position (atoms!!2))
+> newMolecule name atoms = (name,getBasis atoms,atoms) 
+
+Apply a function to the atoms of the molecule.
+
+> doToMolecule :: ([Atom] -> [Atom]) -> Molecule -> Molecule
+> doToMolecule f  (string,axis,atoms) =  newMolecule string (f atoms)
+
+Apply a function to each of the positions of the atoms of the molecule. 
+
+> doToMolecule' :: (Vec -> Vec) -> (Molecule -> Molecule)
+> doToMolecule' f   =  doToMolecule (map (doToAtom f))
+

+ ZEBEDDE/Core/ReadXYZ.hs view

@@ -0,0 +1,61 @@+{-# LANGUAGE TypeSynonymInstances #-}
++{-# LANGUAGE FlexibleInstances #-}
+module ZEBEDDE.Core.ReadXYZ (readXYZ,writeXYZ,showMolecule) where
+import ZEBEDDE.Core.Vector
+import ZEBEDDE.Core.Molecule
+
+type XYZ = (Int,[(String,Vec)])
+
+readXYZ :: String -> IO Molecule
+readXYZ filepath = do
+               dat <- readFile filepath
+               let n = read (firstLine dat) in do
+                return (newMolecule n (map getName (atoms dat)))
+ where
+  getName (a,b) | a == "C" = (a,6,b)
+                | a == "H" = (a,1,b)
+                | otherwise = (a,0,b) -- this is still to be implemented and should decide the atomic weight from the name. 
+  firstLine = takeWhile (not . (=='\n'))
+  rest a = drop ((length (firstLine a))+2) a
+  atoms a =  map getAtom $ filter ((4==).length ) $ map ((split ' ')) (split '\n' (rest a))
+   where
+    getAtom :: [String] -> (String,Vec)
+    getAtom xss = let [a,b,c,d] = take 4 xss in (a,vec (read b,read c,read d)) 
+    split flag string | string == [] = []
+                      | string == [flag] = []
+                      | (head string) == flag = split flag (tail string)
+                      | flag `elem` string = let x = (takeWhile (not.(==flag)) string) in x:(split flag (drop ((length x)+1) string))
+                      | otherwise = [string] 
+
+writeXYZ :: String -> Molecule -> IO ()
+writeXYZ s xyz = writeFile s (showMolecule xyz) 
+
+showMolecule m = let atoms' = (atoms m) in 
+   (show (length atoms')) +++  (foldl (+++) "" (map show'' atoms'))
+    where
+    (+++) a b = a ++ "\n"++ ( b)
+    show'' ((s,_, v)) = (s) ++ show v
+
+instance Show Vec where
+ show v = (show' ((unVec v)::(Double,Double,Double)) )
+  where 
+   show'  ((a,b,c)) = concat $  map (show''.show''') [a,b,c]
+   show'' a = if (head a) == '-' then ("    " ++ (pad a)) else ("     " ++ (pad a))
+   show''' a = show a
+   pad   a | length a == 11 = a 
+           | otherwise = a ++ (take (11-(length a)) (repeat ' '))
+
+
+
+test = do
+ putStr "please run this in the same directory as \"benzene.xyz\"\n"
+ _ <- getLine
+ xyz <- readXYZ "benzene.xyz"
+ writeXYZ "out.xyz" xyz  
+ putStrLn "thank you, that worked."
+
+
+
+
+

+ ZEBEDDE/Core/Vector.lhs view

@@ -0,0 +1,81 @@+
+\section{Vectors}
+\begin{verbatim}
+           This module exports the type Vec.                              
+           A constructor and destructor for this type are exposed.   
+           Vec is made an instance of Eq,Ord,Num and Fractional;         	
+            meaning functions defined over these types 			
+            are made available for Vec. e.g. sum (a::[Vec])		 			
+\end{verbatim}
+
+> {-# LANGUAGE GeneralizedNewtypeDeriving #-} 
+> module ZEBEDDE.Core.Vector 
+> 			    ( Vec(unVec),
+> 			      vec,
+> 			      apply,
+> 			      apply2,
+> 			      magnitude,
+> 			      testVectorHs) where
+
+
+Create vector type.
+
+> newtype Vec  = Vec  {
+>                   unVec :: (Double,Double,Double)
+>                } deriving (Eq)
+> vec = Vec 
+
+\begin{verbatim}
+           We do not export the newtype Vec datatype because we 
+            wish to restrict patern matching 
+            incase its structure changes later on. 
+           Instead we export the constructor vec, whch cannot be pattern matched. 
+\end{verbatim}\\
+
+Apply a function to each element of a vector
+
+> apply f (Vec (a,b,c)) = Vec (f a,f b,f c)
+
+\begin{verbatim}
+Transform  a function taking two doubles and returning a double,
+into       a function taking two vectors and returning a vector
+\end{verbatim}
+
+> apply2 ::  (Double -> Double -> Double) -> (Vec -> Vec -> Vec)
+> apply2 f (Vec (x,y,z)) (Vec (a,b,c)) = Vec (f x a,f y b, f z c)
+
+Define magnetude of vector
+
+> magnitude (Vec (a,b,c)) =  sqrt (a^2+b^2+c^2)
+
+Order vectors by their magnitude
+
+> instance Ord Vec where
+>  (<=) a b = (magnitude a) <= (magnitude b)
+
+Method to add, subtract and multiply vectors, \\
+also extends some standard functions on the integers to vectors.
+
+> instance Num Vec where
+>  fromInteger a = let a' = fromInteger a in vec (a',a',a')
+>                          -- creates a vector whoes
+>                          -- elements are all the given integer, expressed as doubles.
+>  (+) = apply2 (+)
+>  (-) = apply2 (-)
+>  (*) = apply2 (*)        -- so can e.g. multiply a vector by a scalar simply using (*)  
+>  abs = apply abs         -- absolute value of each element of the vector as a vector. 
+>  signum = apply signum   -- signum of each element of the vector as a vector. 
+
+Create fractional vectors of rational numbers.
+
+\begin{verbatim}
+           While we store Vec as a tripple of type Double,
+            here the actual tripple is made to be a rational number 
+            by element-wise division.
+\end{verbatim}
+
+> instance Fractional Vec where
+>  fromRational a = apply (+(fromRational a)) 0 -- can now opperate on rational scalars, e.g. division by a scalar
+>  (/) = apply2 (/)
+
+> testVectorHs = Vec (0.5,0.5,0.5) == 1/2

+ ZEBEDDE/Quaternion.hs view

@@ -0,0 +1,26 @@+module ZEBEDDE.Quaternion where
+
+import ZEBEDDE.Core.Vector
+
+import Data.Vect.Float.Base
+import Data.Vect.Float.Util.Quaternion
+
+
+vec'  :: Vec  -> Vec3 
+vec' v = (\ (a,b,c)-> Vec3 (realToFrac a)(realToFrac b)(realToFrac c) ) $ unVec v
+vec'' :: Vec3 -> Vec
+vec'' (Vec3 a b c) = vec ((realToFrac a),(realToFrac b),(realToFrac c))
+
+rotate :: (Vec,Double) -> Vec -> Vec
+rotate (axis,angle) v = vec'' (actU (rotU'' axis angle) (vec' v))
+
+rotU'' :: Vec -> Double -> U
+rotU'' v d = rotU (vec' v) (realToFrac d)
+
+
+testRotation n = rotate (vec (1,0,0),n) (vec (0,0,1))
+
+
+
+
+

+ ZEBEDDE/Translation.lhs view

@@ -0,0 +1,56 @@+\section{Translation}
+
+Changing the vector position of an atom via translation an rotation abut a point
+
+> module ZEBEDDE.Translation 
+>              ( translate,
+>                com,
+>                comFrame,
+>                rotate,
+>                rotateAboutPoint,
+>                rotateInCom) where
+
+> import ZEBEDDE.Core.Vector
+> import ZEBEDDE.Core.Molecule
+> import qualified ZEBEDDE.Quaternion as Q
+
+Define translation for molecules
+
+> translate :: Vec -> (Molecule -> Molecule)
+> translate v = doToMolecule (map (doToAtom (+v)))
+
+
+> rotate :: Vec -> Double -> (Molecule -> Molecule)
+> rotate v d = doToMolecule' (Q.rotate (v,d))
+
+Rotate about a point
+
+> rotateAboutPoint :: Vec -> Vec -> Double -> (Molecule -> Molecule)
+> rotateAboutPoint center axis theta = t3.t2.t1
+>  where
+>   t1 = translate (0-center)
+>   t2 = rotate axis theta
+>   t3 = translate center
+
+Perform a rotation on a Molecule in its center of mass frame
+
+> centerOfMass :: [Atom] -> Vec
+> centerOfMass atoms = (sum x) / (sum ms)
+>  where
+>     x = zipWith (*) ms rs
+>     ms = map fromIntegral ms'
+>     (_,ms',rs) = unzip3 atoms
+
+> com = centerOfMass.atoms 
+
+> comFrame :: Molecule -> Molecule
+> comFrame m@(_,_,atoms) = translate (0-(centerOfMass atoms)) m
+
+> rotateInCom :: Vec -> Double -> Molecule -> Molecule 
+> rotateInCom axis theta m = rotateAboutPoint (com m) axis theta m
+
+
+
+
+
+