-
Notifications
You must be signed in to change notification settings - Fork 1
/
Fibril.hs
242 lines (215 loc) · 10.6 KB
/
Fibril.hs
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
{-# OPTIONS_GHC -funbox-strict-fields #-}
-- | Allows for making a repeat polymer.
module Fibril( Fibril (..)
, makeFibril
, extractFibril
, glueCartesianChain -- move to
, instantiate
, FibrilModel (..)
, makeFibrilModel
, delimitFragSet
, sampleFibrilMonomer
, sampleFibril
, sampleFibrilModel ) where
import Control.Monad.ST
--import Control.Monad.Primitive
import System.Random (RandomGen, randomR)
import qualified Data.Vector as V
import Data.Maybe(fromMaybe, fromJust) -- DEBUG: fromMaybe
import Data.Ix(inRange)
import Control.Exception(assert)
import Control.DeepSeq(NFData(..))
import Data.Vector.V3 as V3
import Util.Angle(degree2radian)
import Util.Random(normal)
import Debug.Trace(traceShow) -- DEBUG
import Rosetta.Fragments( RFragSet(..)
, RFrag (..) )
import Topo
import FragReplace
import Model
-- * Fibril data structure and its expansion.
-- | Holds all information necessary for constructing a polymer
-- of N monomer repeats, without any linker.
data Fibril = Fibril { monomer :: TorsionTopo
, shift :: !Double
, twist :: !Double
, repeats :: !Int
, monomerLen :: !Int
}
deriving(Show)
instance NFData Fibril where
rnf = rnf . monomer
-- | Makes a Polymer out of a monomer and linker topologies, and a number
-- of monomer repeats.
makeFibril monomerTopo repeats shift twist = assert (repeats > 1) $
Fibril {
monomer = monomerTopo'
, shift = shift
, twist = twist
, repeats = repeats
, monomerLen = monoLen
}
where
monomerTopo' = renumberResiduesT 1 monomerTopo
monoLen = lastResidueId monomerTopo'
-- * PolymerModel as used during modeling.
-- | Polymer model data structure and its cached expansions.
data FibrilModel = FModel { fibril :: Fibril
, tFibril :: TorsionTopo -- with "gapped" connections!
, cFibril :: CartesianTopo
}
-- TODO: here TorsionTopo is of different length than CartesianTopo!
instance Model FibrilModel where
cartesianTopo fm = debugging fm $ cFibril fm
where
debugging fm = id
--debugging fm = traceShow ("Fibril::cartesianTopo", length $ filter ((=="CA") . cAtName) $ backbone $ cFibril fm)
torsionTopo = tFibril -- there are no "cut" nodes in TorsionTopo so far!
instance NFData FibrilModel where
rnf m = rnf (fibril m) `seq` rnf (cFibril m)
-- | Converts Polymer to a PolymerModel, by lazy caching of its expansions
-- as both TorsionTopo and CartesianTopo.
makeFibrilModel f = FModel { fibril = f
, tFibril = monomer f
, cFibril = instantiate f
}
-- | Instantiates a polymer of N identical monomer units, linked by a linker unit.
instantiate :: Fibril -> CartesianTopo
instantiate f = debugging result
where
result = (renumberAtomsC 1 $
renumberResiduesC 1 $
foldr1 glueCartesianChain $
zipWith shiftMonomer [1..repeats f] $
repeat $
computePositions $
monomer f )
shiftMonomer num mono = fmap xform mono
where
-- TODO: add twist!
xform cart@(Cartesian { cPos = pos }) =
cart { cPos = shiftTwist (shift f * fromIntegral num)
(twist f ) pos }
debugging = id
{-debugging = traceShow ("instantiate",
length $ filter ((=="CA") . cAtName) $ backbone result,
repeats f) -}
{-# INLINE shiftTwist #-}
-- | Shift & twist along Z axis.
shiftTwist s t = shiftZ s . twistXY t
-- | Translation along Z axis.
shiftZ s (V3.Vector3 x y z) = V3.Vector3 x y $ z + s
-- | Counterclockwise twist along Z axis (along XY plane.)
twistXY t (V3.Vector3 x y z) = V3.Vector3 (x*cost - y*sint)
(x*sint + y*cost) z
where
sint = sin $ degree2radian t
cost = cos $ degree2radian t
-- | Takes TorsionTopo of a polymer, starting, ending residues of a monomer,
-- and an ending residue of a linker afterwards, and constructs an N-repeat
-- polymer. That means that we cannot pick last monomer, and broken torsion
-- angles suggest we should pick first monomer either.
extractFibril :: Int -> Int -> Int -> TorsionTopo -> Either String Fibril
extractFibril i j n topo = do (_, _, monomerAndTail) <- tagError "Cannot find start of the monomer" $
topo `splitTopoAt` findResId i
(monomer, tail ) <- tagError "Cannot find end of the monomer" $
monomerAndTail `cutTopoAt` findResId (j+1)
return $! makeFibril monomer n 4.8 0.0
-- TODO: extract shift&twist
where
-- TODO: tagError may be generic transformation from Maybe monad to Either monad! (Or `ifMissing`?)
tagError :: String -> Maybe a -> Either String a
tagError msg = maybe (Left msg) Right
-- | Checks that we are at given residue id.
findResId i torsion = tResId torsion == i
-- * Utilities for gluing chains
-- | Replaces OXT of first chain with the second chain.
-- TODO: move it to FragReplacement.
glueCartesianChain :: CartesianTopo -> CartesianTopo -> CartesianTopo
topo1 `glueCartesianChain` topo2 = debugging result
where
Just result = topo1 `changeTopoAt` (\t -> cAtName t == "OXT") $ const topo2
debugging = id
{- debugging = traceShow ("glueCartesianChain",
length $ backbone topo1 ,
length $ backbone topo2 ,
length $ backbone result) -}
-- | Prepare fragment set by removing those that span the boundary between
-- monomer and linker.
delimitFragSet :: Fibril -> RFragSet -> RFragSet
delimitFragSet fibril = RFragSet . V.filter criterion . unRFragSet
where
notOnBoundary :: RFrag -> Bool
notOnBoundary f = endPos f <= monomerLen fibril
criterion :: V.Vector RFrag -> Bool
criterion = notOnBoundary . V.head
-- | Sampling of a polymer with random fragment.
sampleFibrilMonomer fragSet fibril gen = (fibrilXchgMono, gen')
where
fibrilXchgMono = fibril { monomer = xchg $ monomer fibril }
-- TODO: change fromMaybe to fromJust, after cutting fragment set at right place
xchg topo = fromJust $ replaceFragment pos frag topo
((_pos, frag), gen') = randomF fragSet gen
pos = startPos frag
assertions = assert $ pos <= monomerLen fibril
-- | Stationary process along the given mean, with given standard deviation (for sampling shift and twist.)
transition :: (RandomGen r) => Double -- ^ target (base) value
-> Double -- ^ standard deviation of a single sampling
-> Double -- ^ current (starting) value
-> r -- ^ random number generator
-> (Double, r)
transition mean stdev current gen = normal ((mean + current)/2.0) stdev gen
-- NOTE: Standard deviation at starting temperature of 1.0.
baseTwist = 0
baseTwistDev = 30 -- +-30 degrees
baseShift = 4.8
baseShiftDev = 0.8 -- 0.8 Angstroem
-- current * exp (-theta * t) + mean (1-exp(-theta*t) + stdev/sqrt (2
-- | Sampling method to be applied to whole FibrilModel.
--
-- TODO: Sampling of twist and shift should depend on temperature?
sampleFibril :: RandomGen t => Double -- ^ probability of the shift transition.
-> Double -- ^ probability of the twist transition.
-> RFragSet -- ^ fragment set for replacement transition.
-> Fibril -- ^ input Fibril.
-> t -- ^ random number generator
-> (Fibril, t)
sampleFibril shiftProb twistProb fragSet aFibril gen = assertions
(newFibril, gen'')
where
-- Pick type of MC move.
(shiftOrTwist, gen') = randomR (0.0, 1.0) gen
moveFragment = shiftOrTwist >= shiftProb + twistProb
moveShift = shiftOrTwist <= shiftProb
moveTwist = shiftOrTwist > shiftProb && shiftOrTwist <= twistProb
-- Twist sampling
(newTwist, gen''twist) = transition baseTwist baseTwistDev (twist aFibril) gen'
-- Shift sampling
(newShift, gen''shift) = transition baseShift baseShiftDev (shift aFibril) gen'
-- Perform MC move.
(newFibril, gen'') =
if moveShift
then (aFibril { shift = newShift }, gen''shift) -- TODO: implement
else
if moveTwist
then (aFibril { twist = newTwist }, gen''twist) -- TODO: implement
else assert moveFragment $
sampleFibrilMonomer fragSet aFibril gen'
-- Check input parameters
assertions x = assert (shiftProb + twistProb < 1.0) $
assert (shiftProb > 0.0) $
assert (twistProb > 0.0) $
assert (moveShift /= (moveTwist /= moveFragment)) $
assert (moveFragment <= not moveShift ) x
-- | Sampling function for a fibril model, either changes shift&twist parameters,
-- or exchanges a fragment within a monomer.
sampleFibrilModel :: RandomGen t => Double -- ^ Probability of the shift transition.
-> Double -- ^ Probability of the twist transition.
-> RFragSet -- ^ Fragment set for replacement transition.
-> FibrilModel -- ^ Input model.
-> t -- ^ Random number generator
-> (FibrilModel, t)
sampleFibrilModel shiftProb twistProb fragSet fm gen = (makeFibrilModel newFibril, gen')
where
(newFibril, gen') = sampleFibril shiftProb twistProb fragSet (fibril fm) gen