← Back to branch summary

~ubuntu-branches/debian/sid/lammps/sid

« back to all changes in this revision

Viewing changes to tools/moltemplate/examples/CG_biomolecules/metaphase_chromatin_Naumova+Imakaev2013/moltemplate_files/generate_system_lt.py

  • Committer: Package Import Robot
  • Author(s): Anton Gladky
  • Date: 2015-04-29 23:44:49 UTC
  • mfrom: (5.1.3 experimental)
  • Revision ID: package-import@ubuntu.com-20150429234449-mbhy9utku6hp6oq8
Tags: 0~20150313.gitfa668e1-1
Upload into unstable.

Show diffs side-by-side

added added

removed removed

Lines of Context:
tools/moltemplate/examples/CG_biomolecules/metaphase_chromatin_Naumova+Imakaev2013/moltemplate_files/generate_system_lt.py
1
 
#!/usr/bin/env python
2
 
 
3
 
err_msg = """
4
 
Usage:
5
 
 
6
 
   generate_system_lt.py n < monomer_coords.raw > system.lt
7
 
 
8
 
Example:
9
 
 
10
 
   generate_system_lt.py 30118 47 < coords.raw > system.lt
11
 
 
12
 
Explanation:
13
 
     ARGUMENTS:
14
 
         n = total length of the polymer (in monomers)
15
 
         L = the (average) length of each condensin interval (Poisson-
16
 
             distributed)  This is also 1/probability that each monomer
17
 
             is a "condensin monomer".
18
 
 
19
 
   (Note: 30117 ~= 128000/4.25, but using 30118 makes interpolation cleaner,
20
 
          and 47 = 200/4.25. Note that 128000 and 200 are for the 10nm model.
21
 
          See the supplemental section of Naumova et al Science 2013, p 18.)
22
 
 
23
 
"""
24
 
 
25
 
 
26
 
import sys
27
 
import random
28
 
from math import *
29
 
 
30
 
# Parse the argument list:
31
 
if len(sys.argv) <= 2:
32
 
    sys.stderr.write("Error:\n\nTypical Usage:\n\n"+err_msg+"\n")
33
 
    exit(1)
34
 
N=int(sys.argv[1])
35
 
L=float(sys.argv[2])
36
 
if len(sys.argv) > 3:
37
 
    delta_x = float(sys.argv[3])
38
 
else:
39
 
    delta_x = 2.0
40
 
if len(sys.argv) > 4:
41
 
    x_offset = float(sys.argv[4])
42
 
else:
43
 
    x_offset = -((N-1.0)/2) * delta_x
44
 
 
45
 
 
46
 
coords = [[0.0, 0.0, 0.0] for i in range(0,N)]
47
 
lines = sys.stdin.readlines()
48
 
if len(lines) != N:
49
 
    sys.stderr.write("Error: Number of lines in input file ("+str(len(lines))+")\n"
50
 
                     "       does not match first argument ("+str(N)+")\n")
51
 
    exit(1)
52
 
for i in range(0, N):
53
 
    coords[i] = map(float, lines[i].split())
54
 
 
55
 
# Now calculate the box_boundaries:
56
 
box_bounds_min = [0.0, 0.0, 0.0]
57
 
box_bounds_max = [0.0, 0.0, 0.0]
58
 
for i in range(0, N):
59
 
    for d in range(0, 3):
60
 
        if i == 0:
61
 
            box_bounds_min[d] = coords[i][d]
62
 
            box_bounds_max[d] = coords[i][d]
63
 
        else:
64
 
            if coords[i][d] > box_bounds_max[d]:
65
 
                box_bounds_max[d] = coords[i][d]
66
 
            if coords[i][d] < box_bounds_min[d]:
67
 
                box_bounds_min[d] = coords[i][d]
68
 
 
69
 
# Now scale the box boundaries outward by 50%
70
 
box_scale = 1.5
71
 
for d in range(0,3):
72
 
    box_bounds_cen    = 0.5*(box_bounds_max[d] + box_bounds_min[d])
73
 
    box_bounds_width  = box_bounds_max[d] - box_bounds_min[d]
74
 
    box_bounds_min[d] = box_bounds_cen - 0.5*box_scale*box_bounds_width
75
 
    box_bounds_max[d] = box_bounds_cen + 0.5*box_scale*box_bounds_width
76
 
 
77
 
# Now calculate the direction each molecule should be pointing at:
78
 
direction_vects = [[0.0, 0.0, 0.0] for i in range(0,N)]
79
 
for d in range(0, 3):
80
 
     direction_vects[0][d] = coords[1][d] - coords[0][d]
81
 
     direction_vects[N-1][d] = coords[N-1][d] - coords[N-2][d]
82
 
for i in range(1, N-1):
83
 
    for d in range(0, 3):
84
 
         direction_vects[i][d] = coords[i+1][d] - coords[i-1][d]
85
 
 
86
 
# Optional: normalize the direction vectors
87
 
for i in range(1, N-1):
88
 
    direction_len = 0.0
89
 
    for d in range(0, 3):
90
 
        direction_len += (direction_vects[i][d])**2
91
 
    direction_len = sqrt(direction_len)
92
 
    for d in range(0, 3):
93
 
        direction_vects[i][d] /= direction_len
94
 
 
95
 
# Now, begin writing the text for the system.lt file:
96
 
 
97
 
sys.stdout.write(
98
 
"""
99
 
import "monomer.lt"      # <-- defines "Monomer"
100
 
import "condensin.lt"    # <-- defines "CondensinMonomer"
101
 
 
102
 
 
103
 
"""
104
 
)
105
 
 
106
 
 
107
 
 
108
 
# Figure out which monomers are "Monomers" and which monomers are 
109
 
# "CondensinMonomers"
110
 
 
111
 
ic = 0 # count the number of condensins added so far
112
 
condensin_is_here = [False for i in range(0, N)]
113
 
for i in range(0, N):
114
 
    #add_condensin_here = random.random() <  (1.0 / L)
115
 
    add_condensin_here = random.random() <  (1.0 / (L-2.0))
116
 
 
117
 
    # We do not allow condensin at successive sites separated by less than 2
118
 
    # subunits (the "L-2.0" above is to compensate for this)
119
 
    if (((i > 0) and condensin_is_here[i-1]) or
120
 
        ((i > 1) and condensin_is_here[i-2])):
121
 
        add_condensin_here = False
122
 
 
123
 
    if add_condensin_here:
124
 
        condensin_is_here[i] = True
125
 
        ic += 1
126
 
Nc = ic
127
 
 
128
 
 
129
 
ic = 0
130
 
for i in range(0, N):
131
 
    if condensin_is_here[i]:
132
 
        sys.stdout.write("condensins["+str(ic)+"] = new CondensinMonomer.scale(0.5,0.8,0.8).rotvv(1,0,0,")
133
 
        ic+=1
134
 
    else:
135
 
        sys.stdout.write("monomers["+str(i)+"] = new Monomer.scale(0.5,0.8,0.8).rotvv(1,0,0,")
136
 
    sys.stdout.write(str(direction_vects[i][0])+","
137
 
                     +str(direction_vects[i][1])+","
138
 
                     +str(direction_vects[i][2])+
139
 
                     ").move("
140
 
                     +str(coords[i][0])+","
141
 
                     +str(coords[i][1])+","
142
 
                     +str(coords[i][2])+")\n")
143
 
 
144
 
    #if condensin_is_here[i]:
145
 
    #    if i < N-1:
146
 
    #        sys.stdout.write("\n"
147
 
    #                         "#(override the dihedral angle for this monomer)\n"
148
 
    #                         "write(\"Data Dihedrals\") {\n"
149
 
    #                         "  $dihedral:twistor"+str(i+1)+" @dihedral:CondensinMonomer/TWISTOR $atom:monomers["+str(i)+"]/t $atom:monomers["+str(i)+"]/c $atom:monomers["+str(i+1)+"]/c $atom:monomers["+str(i+1)+"]/t\n"
150
 
    #                         "}\n"
151
 
    #                         "\n")
152
 
 
153
 
 
154
 
 
155
 
sys.stdout.write(
156
 
"""
157
 
 
158
 
# ---------------- simulation box -----------------
159
 
 
160
 
# Now define a box big enough to hold a polymer with this (initial) shape
161
 
 
162
 
"""
163
 
)
164
 
 
165
 
 
166
 
sys.stdout.write("write_once(\"Data Boundary\") {\n"
167
 
                 +str(box_bounds_min[0])+"  "+str(box_bounds_max[0])+" xlo xhi\n"
168
 
                 +str(box_bounds_min[1])+"  "+str(box_bounds_max[1])+" ylo yhi\n"
169
 
                 +str(box_bounds_min[2])+"  "+str(box_bounds_max[2])+" zlo zhi\n"
170
 
                 "}\n\n\n")
171
 
 
172
 
 
173
 
sys.stdout.write(
174
 
"""
175
 
# What kind of boundary conditions are we using?
176
 
 
177
 
write_once("In Init") {
178
 
  boundary s s s      # <-- boundary conditions in x y z directions
179
 
  #boundary p p p      # <-- boundary conditions in x y z directions
180
 
}
181
 
# "p" stands for "periodic"
182
 
# "s" stands for "shrink-wrapped" (non-periodic)
183
 
 
184
 
 
185
 
# ---- Bonds ----
186
 
 
187
 
 
188
 
write_once("In Settings") {
189
 
  #  10nm model:
190
 
  #bond_coeff @bond:backbone harmonic 100.0 1.0
191
 
  #  30nm fiber (4.25^(1/3)=1.6198059006387417)
192
 
  bond_coeff @bond:backbone harmonic 100.0 1.6198059006387417
193
 
}
194
 
 
195
 
 
196
 
"""
197
 
)
198
 
 
199
 
 
200
 
sys.stdout.write("write(\"Data Bonds\") {\n")
201
 
 
202
 
# Old bond-loop was simple:
203
 
#for i in range(0, N-1):
204
 
#     sys.stdout.write("  $bond:b"+str(i+1)+" @bond:backbone $atom:monomers["+str(i)+"]/a  $atom:monomers["+str(i+1)+"]/a\n")
205
 
 
206
 
ic = 0
207
 
for i in range(0, N-1):
208
 
    #sys.stderr.write("i="+str(i)+", ic="+str(ic)+", Nc="+str(Nc)+"\n")
209
 
 
210
 
    # Figure out if the first atom in the bond pair 
211
 
    # belongs to a regular Monomer or a CondensinMonomer
212
 
    if condensin_is_here[i]:
213
 
        sys.stdout.write("  $bond:b"+str(i+1)+" @bond:backbone $atom:condensins["+str(ic)+"]/a")
214
 
        ic+=1
215
 
    else:
216
 
        sys.stdout.write("  $bond:b"+str(i+1)+" @bond:backbone $atom:monomers["+str(i)+"]/a")
217
 
 
218
 
    # Do the same thing for the second atom in the bond pair
219
 
    if condensin_is_here[i+1]:
220
 
        assert(ic<Nc)
221
 
        sys.stdout.write("  $atom:condensins["+str(ic)+"]/a\n")
222
 
    else:
223
 
        sys.stdout.write("  $atom:monomers["+str(i+1)+"]/a\n")
224
 
 
225
 
sys.stdout.write("}\n\n\n")
226
 
 
227
 
 
228
 
sys.stdout.write("""
229
 
 
230
 
write_once("Data Angles By Type") {
231
 
  @angle:backbone  @atom:*  @atom:*  @atom:*  @bond:backbone  @bond:backbone
232
 
}
233
 
 
234
 
write_once("In Settings") {
235
 
  # Most parameters here were taken from the supplemental material of 
236
 
  # Naumova et al. Science 2013 (simulations by Maxim Imakaev, see Supp Mat)
237
 
  #angle_coeff @angle:backbone cosine 5.0                 #<-10nm fiber
238
 
  angle_coeff @angle:backbone cosine 1.1764705882352942  #<-30nm fiber
239
 
}
240
 
 
241
 
""")
242
 
 
243
 
 
244
 
sys.stdout.write(
245
 
"""
246
 
 
247
 
# ---- Condensins randomly located on the polymer ----
248
 
 
249
 
# Stage 1:
250
 
# Add bonds between consecutive condensin anchors.
251
 
# Imakaev calls this "stage 1: linear compaction":
252
 
 
253
 
"""
254
 
)
255
 
 
256
 
sys.stdout.write("write(\"Data Bonds\") {\n")
257
 
ic = 0
258
 
for i in range(0, N):
259
 
    if condensin_is_here[i]:
260
 
        if (0 < ic):
261
 
            sys.stdout.write("  $bond:bstage1_"+str(ic-1)+"  @bond:stage1  $atom:condensins["+str(ic-1)+"]/a $atom:condensins["+str(ic)+"]/a\n")
262
 
        ic += 1
263
 
 
264
 
sys.stdout.write("}\n\n")
265
 
 
266
 
 
267
 
 
268
 
sys.stdout.write("""
269
 
 
270
 
# Stage 2:
271
 
# Add additional bonds between all pairs of condensin anchors
272
 
# in a window of |ic-jc| <= 30 anchors (along the chain).
273
 
# In the paper, they call this stage 2 axial compression.
274
 
 
275
 
write("Data Bonds") {
276
 
""")
277
 
 
278
 
jcwindowsize = 30
279
 
for ic in range(0, Nc):
280
 
    #jcmin = max(ic-jcwindowsize, 0)
281
 
    #jcmax = min(ic+jcwindowsize, Nc-1)
282
 
    jcmax = min(ic+jcwindowsize, Nc-1)
283
 
    for jc in range(ic+2, jcmax+1):
284
 
        sys.stdout.write("  $bond:bstage2_"+str(ic)+"_"+str(jc)+"  @bond:stage2  $atom:condensins["+str(ic)+"]/a  $atom:condensins["+str(jc)+"]/a\n")
285
 
 
286
 
sys.stdout.write("}\n")
287
 
 
288
 
 
289
 
sys.stdout.write("""
290
 
 
291
 
write_once("In Settings") {
292
 
  # stage 1 bonds are initially off
293
 
  bond_coeff @bond:stage1 harmonic 0.0 0.5   # <--(we can override this later)"
294
 
  # stage 2 bonds are initially off
295
 
  bond_coeff @bond:stage2 harmonic 0.0 0.0   # <--(we can override this later)"
296
 
}
297
 
 
298
 
""")
299
 
 
300
 
sys.stdout.write("\n\n# "+str(Nc)+" condensin molecules added\n\n")
301
 
sys.stderr.write("\n\n# "+str(Nc)+" condensin molecules added\n\n")