Revision | 073836f160b7daa917b1f49e41d3d7c82beb2360 (tree) |
---|---|
Time | 2007-10-19 19:47:15 |
Author | iselllo |
Commiter | iselllo |
I added the code mytest_3_particle_kind.tcl which simulates the dynamics
of two particles estabilishing a fene bond when they are close.
@@ -0,0 +1,170 @@ | ||
1 | +set n_part 2; set density 0.0007 | |
2 | +set box_l [expr pow($n_part/$density,1./3.)] | |
3 | + | |
4 | +setmd box_l $box_l $box_l $box_l | |
5 | +setmd periodic 1 1 1 | |
6 | + | |
7 | + | |
8 | +inter 0 fene 7.0 4.5 | |
9 | +# I define here a fene potential | |
10 | + | |
11 | +set q 0; set type 0 | |
12 | + | |
13 | +# I now set the particle charge to 0 and see if the code works the same of not | |
14 | + | |
15 | +for {set i 0} { $i < $n_part } {incr i} { | |
16 | +set posx [expr $box_l*[t_random]] | |
17 | +set posy [expr $box_l*[t_random]] | |
18 | +set posz [expr $box_l*[t_random]] | |
19 | +part $i pos $posx $posy $posz q $q type $type | |
20 | +} | |
21 | + | |
22 | + | |
23 | +puts "[part ]" | |
24 | + | |
25 | +setmd time_step 0.01; setmd skin 0.4 | |
26 | +set temp 0.0005; set gamma 0.01 | |
27 | +thermostat langevin $temp $gamma | |
28 | + | |
29 | + | |
30 | + | |
31 | + | |
32 | + | |
33 | +set sig 1.0; set cut [expr 1.12246*$sig] | |
34 | +set eps 1.0; set shift [expr 0.25*$eps] | |
35 | +#inter 0 0 lennard-jones $eps $sig $cut $shift 0 | |
36 | +#inter 0 0 morse 1000. 0.01 1. 1.3 | |
37 | + | |
38 | + | |
39 | +#probably the warm-up is now useless since there is no lennard-jones potential | |
40 | + | |
41 | + | |
42 | + | |
43 | + | |
44 | + | |
45 | + | |
46 | + | |
47 | +set integ_steps 200 | |
48 | + | |
49 | +for {set cap 20} {$cap < 200} {incr cap 20} { | |
50 | +puts "t=[setmd time] E=[analyze energy total]" | |
51 | +inter ljforcecap $cap; integrate $integ_steps | |
52 | +set min [analyze mindist] | |
53 | + | |
54 | +} | |
55 | + | |
56 | +puts "Warmup finished. Minimal distance now $min" | |
57 | +#uncap forces | |
58 | +inter ljforcecap 0 | |
59 | + | |
60 | +puts "end of equilibration" | |
61 | + | |
62 | +set vmd "yes" | |
63 | + | |
64 | +if { $vmd == "yes" } { | |
65 | +# This calls a small tcl script which starts the program # | |
66 | +# VMD and opens a socket connection between ESPResSo and # | |
67 | +# VMD. # | |
68 | + prepare_vmd_connection tutorial 3000 | |
69 | + | |
70 | +# Just wait a moment until VMD has started. # | |
71 | +# The 'exec' command is quite useful since with that you can# | |
72 | +# call any other program from within your simulation script.# | |
73 | + exec sleep 4 | |
74 | + | |
75 | +# The additional command imd steers the socket connection # | |
76 | +# to VMD, e.g. sending the actual coordinates # | |
77 | + imd positions | |
78 | +} | |
79 | + | |
80 | + | |
81 | +#prepare the saving of the results | |
82 | + | |
83 | +#set obs [open "rg_lorenzo.dat" "w"] | |
84 | + | |
85 | + | |
86 | +proc and {a b} {return [expr $a && $b]} | |
87 | + | |
88 | +set activate_fene 1 | |
89 | +set fene_warning 1 | |
90 | +set integ_step 20 | |
91 | + | |
92 | +for {set i 0} { $i < 2000 } { incr i} { | |
93 | +set temp [expr [analyze energy kinetic]/(1.5*$n_part)] | |
94 | +#puts "t=[setmd time] E=[analyze energy total], T=$temp" | |
95 | +integrate $integ_steps | |
96 | + | |
97 | +# puts "save configuration" | |
98 | + | |
99 | + | |
100 | +#puts "[analyze nbhood 1 2.]" | |
101 | + | |
102 | + | |
103 | +set acti [analyze nbhood 1 4.5] | |
104 | + | |
105 | + | |
106 | +#puts "acti is, $acti" | |
107 | +set myfene [llength $acti] | |
108 | +if {$myfene > 1 && $fene_warning == 1} { | |
109 | +puts "the particles are within the FENE range" | |
110 | +set fene_warning 0 | |
111 | +} | |
112 | + | |
113 | +set neighbor [analyze nbhood 1 2.] | |
114 | +set mylen [llength $neighbor] | |
115 | +#puts "mylen is, $mylen" | |
116 | + | |
117 | +if { $mylen == 2 && $activate_fene == 1 } { | |
118 | +part 0 bond 0 1 | |
119 | +puts "activated fene bond!" | |
120 | +set activate_fene 0 | |
121 | +} | |
122 | +#part 0 bond 0 1 | |
123 | +#set activate_fene 0 | |
124 | + | |
125 | + | |
126 | + | |
127 | + | |
128 | +set f [open "config_$i" "w"] | |
129 | +blockfile $f write tclvariable {box_l density} | |
130 | +blockfile $f write variable box_l | |
131 | +blockfile $f write particles {id pos type} | |
132 | +close $f | |
133 | +# if you have turned on the vmd option you can now | |
134 | + # follow what your simulation is doing | |
135 | + if { $vmd == "yes" } { imd positions } | |
136 | + | |
137 | +#Now I calculate the radius of gyration | |
138 | +# set rg [lindex [analyze rg] 0] | |
139 | +#puts $obs "[setmd time] $rg" | |
140 | + | |
141 | +} | |
142 | + | |
143 | + | |
144 | + | |
145 | + | |
146 | + | |
147 | +#close $obs | |
148 | +puts "end of integration" | |
149 | + | |
150 | +#plotObs "rg.dat" { 1:2 } labels { "time" "rg" } out "rg" | |
151 | +#exec gv rg.ps | |
152 | + | |
153 | +set f [open "config_1" "r"] | |
154 | +while { [blockfile $f read auto] != "eof" } {} | |
155 | +close $f | |
156 | + | |
157 | +puts "ok reading the block file" | |
158 | + | |
159 | +set rdf [analyze rdf 0 0 0.9 [expr $box_l/2] 100] | |
160 | +set rlist "" | |
161 | +set rdflist "" | |
162 | +foreach value [lindex $rdf 1] { | |
163 | +lappend rlist [lindex $value 0] | |
164 | +lappend rdflist [lindex $value 1] | |
165 | +} | |
166 | + | |
167 | + | |
168 | + | |
169 | + | |
170 | +puts "So far so good" | |
\ No newline at end of file |