提交: 8a76795b249604fea20655be7267d81550c8d7cc - myprojects-hg-reborn (hg)

Commit MetaInfo

• 修訂: 8a76795b249604fea20655be7267d81550c8d7cc (tree)
• 時間: 2008-09-06 00:07:55
• 作者: iselllo
• Commiter: iselllo

Log Message

I added a code for studying aggregate diffusion using the DPD
thermostat.

Change Summary

• add: TCL-codes/restart_aggregates_dpd.tcl

差異

diff -r f42379b7664f -r 8a76795b2496 TCL-codes/restart_aggregates_dpd.tcl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/TCL-codes/restart_aggregates_dpd.tcl	Fri Sep 05 15:07:55 2008 +0000
@@ -0,0 +1,365 @@
+#I now read a single configuration already containing many clusters
+
+# number of monomers in each cluster
+
+set n_part 50
+
+#number of clusters
+
+set type_of_part 200
+
+#total number of monomers
+
+set tot_part [expr $n_part*$type_of_part]
+
+puts "tot_part is, $tot_part"
+
+#Now the situation is different: the box size is chosen and density is not relevant.
+
+
+set box_l 10000.
+
+# set density 0.03
+
+# puts "the density is, $density"
+
+# set volume [expr pow($box_l,3.)]
+
+
+
+
+set magnification 1.
+
+setmd box_l [expr $magnification*$box_l] [expr $magnification*$box_l] [expr $magnification*$box_l]
+setmd periodic 1 1 1
+
+puts "the box side is, $box_l"
+
+
+set q 0
+
+
+
+
+
+# set j 14
+
+
+# Now I am going to read a simple text file, NOT a block file
+#set filename "clusters_ini_$j"
+set filename "copy_of_cluster50_200" 
+#set filename "initial_state_folded" 
+
+set fp [open $filename "r"]
+set data [split [read $fp [file size $filename]] "\n"]
+close $fp
+
+
+
+
+
+ for {set i 0} { $i < $tot_part } {incr i} {
+set p [lindex $data $i]
+#puts "p is, $p" 
+#set p2 [expr $box_l*$p]
+#puts "p2 is, $p2"
+#set pos_part $p
+#puts "pos_part[1] is, $pos_part(1)"
+
+set type [expr $i/$n_part]
+
+puts "type is, $type"
+ 
+
+set posx [lindex $p 0 ]
+set posy [lindex $p 1 ]
+set posz [lindex $p 2 ]
+
+
+# set posx [expr $posx*$box_l]
+# set posy [expr $posy*$box_l]
+# set posz [expr $posz*$box_l]
+
+#puts "posz is, $posz"
+
+part $i pos $posx $posy $posz q $q type $type
+
+
+ }
+
+
+
+
+#parameters of the simulation
+
+
+set tot_time 400.
+
+
+set my_step 0.00125
+
+
+set N_step [expr $tot_time/$my_step]
+
+set N_step [expr round($N_step)]
+
+
+set integ_steps 1
+
+# number of configurations I want to save
+
+set N_save 4000
+#I save a configuration every time steps
+
+set every [expr $N_step/$N_save]
+
+puts "every is, $every"
+
+
+
+puts "the total number of time steps is, $N_step"
+
+
+setmd time_step $my_step; setmd skin 0.4
+set temp 0.5 ; set gamma  1. 
+
+
+set cutoff 1.06
+
+set WF 1
+
+
+set TGAMMA 0.
+
+set TCUTOFF 0.
+
+set TWF 1
+
+
+# thermostat langevin $temp $gamma
+# Now I am going to introduce a different thermostat. See email with Christoph Junghans
+
+thermostat inter_dpd $temp
+
+
+
+#setting up the interaction
+
+
+
+set test_pos [part 0 print pos]
+
+puts "the position of the zero-th particle is, $test_pos"
+
+# Now I can set up the interaction
+
+
+for {set i 0} { $i <= $type_of_part } { incr i} {
+
+inter $i $i tabulated tabulated_interaction.dat
+
+# Now I also add another interaction potential, which come directly from the DPD thermostat 
+# inter X Y inter_dpd GAMMA CUTOFF WF TGAMMA TCUTOFF TWF
+
+inter $i $i inter_dpd $gamma $cutoff $WF $TGAMMA $TCUTOFF $TWF
+
+
+
+}
+
+
+# I am not sure if the following part is really needed; a non-defined interaction between particles of kind i and j should
+# be equivalent to setting it to zero.
+
+
+# #Now I explicitely set the dpd interaction as zero for monomers of two different kinds
+
+
+# for {set i 0} { $i <= $type_of_part } { incr i} {
+
+# for {set j 0} { $j <= $type_of_part } { incr j} {
+
+#     if {$i != $j} {
+
+# # inter $i $i inter_dpd gamma cutoff WF TGAMMA TCUTOFF TWF
+
+# inter $i $j inter_dpd 0. 0. 1 0. 0. 1
+
+
+# }
+
+# }
+
+# }
+
+
+
+# }
+
+
+
+
+
+
+
+
+
+set vmd "no"
+
+
+puts "OK here after inter_dpd"
+
+
+if { $vmd == "yes" } {
+# This calls a small tcl script which starts the program    #
+# VMD and opens a socket connection between ESPResSo and    #
+# VMD.                                                      #
+    prepare_vmd_connection tutorial 3000
+
+# Just wait a moment until VMD has started.                 #
+# The 'exec' command is quite useful since with that you can#
+# call any other program from within your simulation script.#
+    exec sleep 4
+
+
+#puts $obs4 "[setmd time] [part 0 print pos]"
+
+}
+
+
+
+#prepare the saving of the results
+
+set count 0
+set label 0
+
+
+#set obs [open "temp.dat" "w"]
+#set obs2 [open "tot_energy.dat" "w"]
+#set obs3 [open "aggregation.dat" "w"]
+#set obs3 [open "rgyr.dat" "w"] 
+
+#set obs4 [open "part_pos.dat" "w"] 
+set obs5 [open "time.dat" "w"] 
+
+
+#initialization random numbers
+
+expr srand(200)
+set l ""; for {set i 0} {$i < [setmd n_nodes]} {incr i} {
+ lappend l [expr int(32768*rand())]
+}
+eval t_random seed $l
+
+
+
+
+#puts "the simulation time now is, [setmd time]"
+
+# for {set i 0} { $i < 3000 } { incr i} {
+for {set i 0} { $i < $N_step } { incr i $integ_steps} {
+
+
+set my_rem [expr $count%$every]
+
+#puts "my_rem is, $my_rem"
+
+
+if {$my_rem== 0 } {
+#set f [open "config_vel_$i" "w"]
+set f [open "config_vel_$label" "w"]
+incr label
+
+#blockfile $f write particles "id folded_position v"
+
+#I use the variant above if I want to save the folded
+#positions, but I can do that myself and it is not interesting
+
+blockfile $f write particles "id pos v"
+
+close $f
+
+
+#set temp [expr [analyze energy kinetic]/(1.5*$n_part)]
+puts "t=[setmd time] "
+
+
+#set energy [lindex [analyze energy total] 0]
+#puts $obs "[setmd time] $temp"
+
+#puts $obs2 "[setmd time] $energy"
+
+
+
+
+#Now I calculate the radius of gyration
+# set rg [lindex [analyze rg 0 1 $n_part] 0]
+#puts "the radius of gyration is, $rg"
+#puts $obs3 "[setmd time] $rg"
+
+puts $obs5 "[setmd time]"
+
+
+# if you have turned on the vmd option you can now
+    # follow what your simulation is doing
+
+# The additional command imd steers the socket connection   #
+# to VMD, e.g. sending the actual coordinates               #
+#    imd positions 
+#pbc wrap -sel "index 13"
+
+
+    if { $vmd == "yes" } { imd positions }
+
+
+
+
+
+
+
+
+
+
+}
+
+incr count
+
+integrate $integ_steps
+
+
+
+#puts $obs4 "[setmd time] [part 0 print pos]"
+
+}
+
+#close $obs
+#close $obs
+#close $obs2
+#close $obs3
+close $obs5
+puts "end of integration" 
+
+puts "the simulation time now is, [setmd time]"
+
+#plotObs "rg.dat" { 1:2 } labels { "time" "rg" } out "rg"
+#exec gv rg.ps
+
+#set f [open "config_1" "r"]
+#while { [blockfile $f read auto] != "eof" } {}
+#close $f
+
+#puts "ok reading the block file"
+
+#set rdf [analyze rdf 0 0 0.9 [expr $box_l/2] 100]
+#set rlist ""
+#set rdflist ""
+#foreach value [lindex $rdf 1] {
+#lappend rlist [lindex $value 0]
+#lappend rdflist [lindex $value 1]
+#}
+
+
+
+
+
+puts "So far so good"
\ No newline at end of file