修訂 | 9e91fa3006a3c1cb4a7a8989966bac6d422a68fe (tree) |
---|---|
時間 | 2007-11-19 21:29:56 |
作者 | iselllo |
Commiter | iselllo |
I modified mytest_1_particle_kind.tcl (cleaned up a bit the code, but
not 100% sure of what I did...in case revert to the previous version to
find out what happened) and added test.tcl, which evolves and saves some
statistics for a set of Brownian particles.
@@ -1,4 +1,4 @@ | ||
1 | -set n_part 50; set density 0.02 | |
1 | +set n_part 50; set density 0.06 | |
2 | 2 | set box_l [expr pow($n_part/$density,1./3.)] |
3 | 3 | |
4 | 4 |
@@ -21,7 +21,7 @@ | ||
21 | 21 | } |
22 | 22 | |
23 | 23 | setmd time_step 0.01; setmd skin 0.4 |
24 | -set temp 0.01; set gamma 0.25 | |
24 | +set temp 0.02; set gamma 0.25 | |
25 | 25 | thermostat langevin $temp $gamma |
26 | 26 | |
27 | 27 |
@@ -31,7 +31,7 @@ | ||
31 | 31 | #set sig 1.0; set cut [expr 1.12246*$sig] |
32 | 32 | #set eps 1.0; set shift [expr 0.25*$eps] |
33 | 33 | |
34 | -set sig 1.0 | |
34 | +set sig 0.1 | |
35 | 35 | set cut 2. |
36 | 36 | set eps 10. |
37 | 37 | set shift 0. |
@@ -0,0 +1,214 @@ | ||
1 | +set n_part 2000; #set density 0.00002 | |
2 | +#set box_l [expr pow($n_part/$density,1./3.)] | |
3 | + | |
4 | +set box_l 2000. | |
5 | +set volume [expr pow($box_l,3.)] | |
6 | + | |
7 | +set density [expr $n_part/$volume ] | |
8 | + | |
9 | +puts "the density is, $density" | |
10 | + | |
11 | + | |
12 | +set magnification 1. | |
13 | + | |
14 | +setmd box_l [expr $magnification*$box_l] [expr $magnification*$box_l] [expr $magnification*$box_l] | |
15 | +setmd periodic 1 1 1 | |
16 | + | |
17 | +puts "the box side is, $box_l" | |
18 | + | |
19 | + | |
20 | +set q 0; set type 0 | |
21 | + | |
22 | +# I now set the particle charge to 0 and see if the code works the same of not | |
23 | + | |
24 | +for {set i 0} { $i < $n_part } {incr i} { | |
25 | +#set posx [expr $box_l*[t_random]] | |
26 | +#set posy [expr $box_l*[t_random]] | |
27 | +#set posz [expr $box_l*[t_random]] | |
28 | + | |
29 | +set posx [expr $box_l/2.] | |
30 | +set posy [expr $box_l/2.] | |
31 | +set posz [expr $box_l/2.] | |
32 | + | |
33 | + | |
34 | +part $i pos $posx $posy $posz q $q type $type | |
35 | +} | |
36 | + | |
37 | + | |
38 | + | |
39 | + | |
40 | +set tot_time 200. | |
41 | + | |
42 | + | |
43 | +set my_step 1. | |
44 | + | |
45 | + | |
46 | +set N_step [expr $tot_time/$my_step] | |
47 | + | |
48 | +set N_step [expr round($N_step)] | |
49 | + | |
50 | + | |
51 | +set integ_steps 1 | |
52 | + | |
53 | + | |
54 | + | |
55 | + | |
56 | +puts "the total number of time steps is, $N_step" | |
57 | + | |
58 | + | |
59 | +setmd time_step $my_step; setmd skin 0.4 | |
60 | +set temp 0.1 ; set gamma 0.1 | |
61 | +thermostat langevin $temp $gamma | |
62 | + | |
63 | + | |
64 | + | |
65 | + | |
66 | + | |
67 | +set sig 1.0 | |
68 | +set cut 20. | |
69 | +set eps 0. | |
70 | +set shift 0. | |
71 | + | |
72 | +#inter 0 0 lennard-jones $eps $sig $cut $shift 0 | |
73 | + | |
74 | + | |
75 | +set tau_increase 10. | |
76 | + | |
77 | +set delta_t [expr $integ_steps*$my_step] | |
78 | + | |
79 | +puts "delta_t is, $delta_t" | |
80 | + | |
81 | +set cap_ini 20. | |
82 | +set cap_fin 200. | |
83 | +set delta_cap [expr $cap_fin-$cap_ini] | |
84 | +set lin_coeff [expr $delta_cap/$tau_increase] | |
85 | + | |
86 | +set n_iter_cap [expr $tau_increase/$delta_t] | |
87 | + | |
88 | +set n_iter_cap [expr round($n_iter_cap)] | |
89 | + | |
90 | +puts "the number of iterations while ramping the potential is, $n_iter_cap" | |
91 | + | |
92 | + | |
93 | +puts "delta_cap is, $delta_cap" | |
94 | + | |
95 | +set cap $cap_ini | |
96 | + | |
97 | + | |
98 | + | |
99 | + | |
100 | +#for {set i 0} { $i <= $n_iter_cap } { incr i} { | |
101 | +#puts "t=[setmd time] E=[analyze energy total]" | |
102 | +#inter ljforcecap $cap; integrate $integ_steps | |
103 | + | |
104 | +#set cap [expr $cap + $lin_coeff*$delta_t ] | |
105 | + | |
106 | +#set min [analyze mindist] | |
107 | + | |
108 | +#} | |
109 | + | |
110 | +#puts "Warmup finished. Minimal distance now $min" | |
111 | +#uncap forces | |
112 | +inter ljforcecap 0 | |
113 | + | |
114 | +puts "end of equilibration" | |
115 | + | |
116 | +set vmd "no" | |
117 | + | |
118 | +if { $vmd == "yes" } { | |
119 | +# This calls a small tcl script which starts the program # | |
120 | +# VMD and opens a socket connection between ESPResSo and # | |
121 | +# VMD. # | |
122 | + prepare_vmd_connection tutorial 3000 | |
123 | + | |
124 | +# Just wait a moment until VMD has started. # | |
125 | +# The 'exec' command is quite useful since with that you can# | |
126 | +# call any other program from within your simulation script.# | |
127 | + exec sleep 4 | |
128 | + | |
129 | +# The additional command imd steers the socket connection # | |
130 | +# to VMD, e.g. sending the actual coordinates # | |
131 | + imd positions | |
132 | +} | |
133 | + | |
134 | + | |
135 | +#prepare the saving of the results | |
136 | + | |
137 | + | |
138 | + | |
139 | + | |
140 | +set obs [open "temp.dat" "w"] | |
141 | +set obs2 [open "tot_energy.dat" "w"] | |
142 | +#set obs3 [open "aggregation.dat" "w"] | |
143 | +set obs3 [open "rgyr.dat" "w"] | |
144 | + | |
145 | +set obs4 [open "part_pos.dat" "w"] | |
146 | + | |
147 | + | |
148 | +set integ_step 20 | |
149 | + | |
150 | + | |
151 | +#puts "the simulation time now is, [setmd time]" | |
152 | + | |
153 | +# for {set i 0} { $i < 3000 } { incr i} { | |
154 | +for {set i 0} { $i < $N_step } { incr i $integ_steps} { | |
155 | + | |
156 | +set temp [expr [analyze energy kinetic]/(1.5*$n_part)] | |
157 | +puts "t=[setmd time] " | |
158 | + | |
159 | + | |
160 | +set energy [lindex [analyze energy total] 0] | |
161 | +puts $obs "[setmd time] $temp" | |
162 | + | |
163 | +puts $obs2 "[setmd time] $energy" | |
164 | + | |
165 | + | |
166 | + | |
167 | +integrate $integ_steps | |
168 | + | |
169 | +#puts "save configuration, [part ]" | |
170 | +part 0 print pos | |
171 | +set f [open "config_$i" "w"] | |
172 | +blockfile $f write tclvariable {box_l density} | |
173 | +blockfile $f write variable box_l | |
174 | +blockfile $f write particles {pos} | |
175 | +close $f | |
176 | +# if you have turned on the vmd option you can now | |
177 | + # follow what your simulation is doing | |
178 | + if { $vmd == "yes" } { imd positions } | |
179 | + | |
180 | +#Now I calculate the radius of gyration | |
181 | + set rg [lindex [analyze rg] 0] | |
182 | +puts $obs3 "[setmd time] $rg" | |
183 | + | |
184 | +puts $obs4 "[setmd time] [part 0 print pos]" | |
185 | + | |
186 | +} | |
187 | + | |
188 | +#close $obs | |
189 | + | |
190 | + | |
191 | +puts "end of integration" | |
192 | + | |
193 | + | |
194 | +puts "the simulation time now is, [setmd time]" | |
195 | + | |
196 | +#plotObs "rg.dat" { 1:2 } labels { "time" "rg" } out "rg" | |
197 | +#exec gv rg.ps | |
198 | + | |
199 | +#set f [open "config_1" "r"] | |
200 | +#while { [blockfile $f read auto] != "eof" } {} | |
201 | +#close $f | |
202 | + | |
203 | +#puts "ok reading the block file" | |
204 | + | |
205 | +#set rdf [analyze rdf 0 0 0.9 [expr $box_l/2] 100] | |
206 | +#set rlist "" | |
207 | +#set rdflist "" | |
208 | +#foreach value [lindex $rdf 1] { | |
209 | +#lappend rlist [lindex $value 0] | |
210 | +#lappend rdflist [lindex $value 1] | |
211 | +#} | |
212 | + | |
213 | + | |
214 | +puts "So far so good" |