修訂 | 7c096a566affaa8b36a889a9b0ccd33b3cb0f22d (tree) |
---|---|
時間 | 2008-11-13 00:15:35 |
作者 | iselllo |
Commiter | iselllo |
I added a code which needs the modified version of Espresso, as it also reads the eta0 and eta1 factors via a
user-supplied array.
@@ -0,0 +1,327 @@ | ||
1 | +#I now read a single configuration already containing many clusters | |
2 | + | |
3 | +# number of monomers in each cluster | |
4 | + | |
5 | +set n_part 10 | |
6 | + | |
7 | +#number of clusters | |
8 | + | |
9 | +set type_of_part 200 | |
10 | + | |
11 | +#total number of monomers | |
12 | + | |
13 | +set tot_part [expr $n_part*$type_of_part] | |
14 | + | |
15 | +puts "tot_part is, $tot_part" | |
16 | + | |
17 | +#Now the situation is different: the box size is chosen and density is not relevant. | |
18 | + | |
19 | + | |
20 | +set box_l 10000. | |
21 | + | |
22 | +# set density 0.03 | |
23 | + | |
24 | +# puts "the density is, $density" | |
25 | + | |
26 | +# set volume [expr pow($box_l,3.)] | |
27 | + | |
28 | + | |
29 | + | |
30 | + | |
31 | +set magnification 1. | |
32 | + | |
33 | +setmd box_l [expr $magnification*$box_l] [expr $magnification*$box_l] [expr $magnification*$box_l] | |
34 | +setmd periodic 1 1 1 | |
35 | + | |
36 | +puts "the box side is, $box_l" | |
37 | + | |
38 | + | |
39 | +set q 0 | |
40 | + | |
41 | + | |
42 | + | |
43 | + | |
44 | + | |
45 | +# set j 14 | |
46 | + | |
47 | + | |
48 | +# Now I am going to read a simple text file, NOT a block file | |
49 | + | |
50 | +set filename "copy_of_cluster10_200" | |
51 | + | |
52 | + | |
53 | +set fp [open $filename "r"] | |
54 | +set data [split [read $fp [file size $filename]] "\n"] | |
55 | +close $fp | |
56 | + | |
57 | + | |
58 | + | |
59 | +#Now I read the eta's for the modified Langevin thermostat | |
60 | + | |
61 | +set filename "eta_0_200" | |
62 | + | |
63 | + | |
64 | +set fp [open $filename "r"] | |
65 | +set eta0 [split [read $fp [file size $filename]] "\n"] | |
66 | +close $fp | |
67 | + | |
68 | + | |
69 | +set filename "eta_1_200" | |
70 | + | |
71 | + | |
72 | +set fp [open $filename "r"] | |
73 | +set eta1 [split [read $fp [file size $filename]] "\n"] | |
74 | +close $fp | |
75 | + | |
76 | + | |
77 | + | |
78 | + | |
79 | + | |
80 | + | |
81 | + | |
82 | + for {set i 0} { $i < $tot_part } {incr i} { | |
83 | +set p [lindex $data $i] | |
84 | + | |
85 | +set eta0_temp [lindex $eta0 $i] | |
86 | + | |
87 | +set eta1_temp [lindex $eta1 $i] | |
88 | + | |
89 | + | |
90 | + | |
91 | + | |
92 | + | |
93 | +set type [expr $i/$n_part] | |
94 | + | |
95 | +puts "type is, $type" | |
96 | + | |
97 | + | |
98 | +set posx [lindex $p 0 ] | |
99 | +set posy [lindex $p 1 ] | |
100 | +set posz [lindex $p 2 ] | |
101 | + | |
102 | + | |
103 | +# set posx [expr $posx*$box_l] | |
104 | +# set posy [expr $posy*$box_l] | |
105 | +# set posz [expr $posz*$box_l] | |
106 | + | |
107 | +#puts "posz is, $posz" | |
108 | + | |
109 | +part $i pos $posx $posy $posz q $q type $type eta $eta0_temp $eta1_temp | |
110 | + | |
111 | + | |
112 | + } | |
113 | + | |
114 | + | |
115 | + | |
116 | +set test_pos [part 0 print pos] | |
117 | + | |
118 | +puts "the position of the zero-th particle is, $test_pos" | |
119 | + | |
120 | +set test_eta [part 0 print eta] | |
121 | + | |
122 | +puts "eta for the zero-th particle is, $test_eta" | |
123 | + | |
124 | + | |
125 | +# Now I can set up the interaction | |
126 | + | |
127 | + | |
128 | +for {set i 0} { $i <= $type_of_part } { incr i} { | |
129 | + | |
130 | +inter $i $i tabulated tabulated_interaction.dat | |
131 | + | |
132 | +} | |
133 | + | |
134 | +#parameters of the simulation | |
135 | + | |
136 | + | |
137 | +set tot_time 100. | |
138 | + | |
139 | + | |
140 | +set my_step 0.00125 | |
141 | + | |
142 | + | |
143 | +set N_step [expr $tot_time/$my_step] | |
144 | + | |
145 | +set N_step [expr round($N_step)] | |
146 | + | |
147 | + | |
148 | +set integ_steps 1 | |
149 | + | |
150 | +# number of configurations I want to save | |
151 | + | |
152 | +set N_save 100 | |
153 | +#I save a configuration every time steps | |
154 | + | |
155 | +set every [expr $N_step/$N_save] | |
156 | + | |
157 | +puts "every is, $every" | |
158 | + | |
159 | + | |
160 | + | |
161 | +puts "the total number of time steps is, $N_step" | |
162 | + | |
163 | + | |
164 | +setmd time_step $my_step; setmd skin 0.4 | |
165 | +set temp 0.5 ; set gamma 1. | |
166 | +thermostat langevin $temp $gamma | |
167 | + | |
168 | + | |
169 | + | |
170 | + | |
171 | +set vmd "no" | |
172 | + | |
173 | + | |
174 | + | |
175 | +if { $vmd == "yes" } { | |
176 | +# This calls a small tcl script which starts the program # | |
177 | +# VMD and opens a socket connection between ESPResSo and # | |
178 | +# VMD. # | |
179 | + prepare_vmd_connection tutorial 3000 | |
180 | + | |
181 | +# Just wait a moment until VMD has started. # | |
182 | +# The 'exec' command is quite useful since with that you can# | |
183 | +# call any other program from within your simulation script.# | |
184 | + exec sleep 4 | |
185 | + | |
186 | + | |
187 | +#puts $obs4 "[setmd time] [part 0 print pos]" | |
188 | + | |
189 | +} | |
190 | + | |
191 | + | |
192 | + | |
193 | +#prepare the saving of the results | |
194 | + | |
195 | +set count 0 | |
196 | +set label 0 | |
197 | + | |
198 | + | |
199 | +#set obs [open "temp.dat" "w"] | |
200 | +#set obs2 [open "tot_energy.dat" "w"] | |
201 | +#set obs3 [open "aggregation.dat" "w"] | |
202 | +#set obs3 [open "rgyr.dat" "w"] | |
203 | + | |
204 | +#set obs4 [open "part_pos.dat" "w"] | |
205 | +set obs5 [open "time.dat" "w"] | |
206 | + | |
207 | + | |
208 | +#initialization random numbers | |
209 | + | |
210 | +expr srand(100) | |
211 | +set l ""; for {set i 0} {$i < [setmd n_nodes]} {incr i} { | |
212 | + lappend l [expr int(32768*rand())] | |
213 | +} | |
214 | +eval t_random seed $l | |
215 | + | |
216 | + | |
217 | + | |
218 | + | |
219 | +#puts "the simulation time now is, [setmd time]" | |
220 | + | |
221 | +# for {set i 0} { $i < 3000 } { incr i} { | |
222 | +for {set i 0} { $i < $N_step } { incr i $integ_steps} { | |
223 | + | |
224 | + | |
225 | +set my_rem [expr $count%$every] | |
226 | + | |
227 | +#puts "my_rem is, $my_rem" | |
228 | + | |
229 | + | |
230 | +if {$my_rem== 0 } { | |
231 | +#set f [open "config_vel_$i" "w"] | |
232 | +set f [open "config_vel_$label" "w"] | |
233 | +incr label | |
234 | + | |
235 | +#blockfile $f write particles "id folded_position v" | |
236 | + | |
237 | +#I use the variant above if I want to save the folded | |
238 | +#positions, but I can do that myself and it is not interesting | |
239 | + | |
240 | +blockfile $f write particles "id pos v" | |
241 | + | |
242 | +close $f | |
243 | + | |
244 | + | |
245 | +#set temp [expr [analyze energy kinetic]/(1.5*$n_part)] | |
246 | +puts "t=[setmd time] " | |
247 | + | |
248 | + | |
249 | +#set energy [lindex [analyze energy total] 0] | |
250 | +#puts $obs "[setmd time] $temp" | |
251 | + | |
252 | +#puts $obs2 "[setmd time] $energy" | |
253 | + | |
254 | + | |
255 | + | |
256 | + | |
257 | +#Now I calculate the radius of gyration | |
258 | +# set rg [lindex [analyze rg 0 1 $n_part] 0] | |
259 | +#puts "the radius of gyration is, $rg" | |
260 | +#puts $obs3 "[setmd time] $rg" | |
261 | + | |
262 | +puts $obs5 "[setmd time]" | |
263 | + | |
264 | + | |
265 | +# if you have turned on the vmd option you can now | |
266 | + # follow what your simulation is doing | |
267 | + | |
268 | +# The additional command imd steers the socket connection # | |
269 | +# to VMD, e.g. sending the actual coordinates # | |
270 | +# imd positions | |
271 | +#pbc wrap -sel "index 13" | |
272 | + | |
273 | + | |
274 | + if { $vmd == "yes" } { imd positions } | |
275 | + | |
276 | + | |
277 | + | |
278 | + | |
279 | + | |
280 | + | |
281 | + | |
282 | + | |
283 | + | |
284 | + | |
285 | +} | |
286 | + | |
287 | +incr count | |
288 | + | |
289 | +integrate $integ_steps | |
290 | + | |
291 | + | |
292 | + | |
293 | +#puts $obs4 "[setmd time] [part 0 print pos]" | |
294 | + | |
295 | +} | |
296 | + | |
297 | +#close $obs | |
298 | +#close $obs | |
299 | +#close $obs2 | |
300 | +#close $obs3 | |
301 | +close $obs5 | |
302 | +puts "end of integration" | |
303 | + | |
304 | +puts "the simulation time now is, [setmd time]" | |
305 | + | |
306 | +#plotObs "rg.dat" { 1:2 } labels { "time" "rg" } out "rg" | |
307 | +#exec gv rg.ps | |
308 | + | |
309 | +#set f [open "config_1" "r"] | |
310 | +#while { [blockfile $f read auto] != "eof" } {} | |
311 | +#close $f | |
312 | + | |
313 | +#puts "ok reading the block file" | |
314 | + | |
315 | +#set rdf [analyze rdf 0 0 0.9 [expr $box_l/2] 100] | |
316 | +#set rlist "" | |
317 | +#set rdflist "" | |
318 | +#foreach value [lindex $rdf 1] { | |
319 | +#lappend rlist [lindex $value 0] | |
320 | +#lappend rdflist [lindex $value 1] | |
321 | +#} | |
322 | + | |
323 | + | |
324 | + | |
325 | + | |
326 | + | |
327 | +puts "So far so good" | |
\ No newline at end of file |