myprojects-hg-reborn

Fork

set n_part 3000; #set density 0.00002
#set box_l [expr pow($n_part/$density,1./3.)]

set box_l 10.
set volume [expr pow($box_l,3.)]

set density [expr $n_part/$volume ]

puts "the density is, $density"


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 type 0

# I now set the particle charge to 0 and see if the code works the same of not

for {set i 0} { $i < $n_part } {incr i} {
#set posx [expr $box_l*[t_random]]
#set posy [expr $box_l*[t_random]]
#set posz [expr $box_l*[t_random]]

#set posx [expr $box_l/2.]
#set posy [expr $box_l/2.]
#set posz [expr $box_l/2.]


set posx 0.
set posy 0.
set posz 0.


part $i pos $posx $posy $posz q $q type $type
}




set tot_time 40.


set my_step 1.


set N_step [expr $tot_time/$my_step]

set N_step [expr round($N_step)]


set integ_steps 1




puts "the total number of time steps is, $N_step"


setmd time_step $my_step; setmd skin 0.4
set temp 0.1 ; set gamma  0.1
thermostat langevin $temp $gamma





set sig 1.0
set cut 20.
set eps 0.
set shift 0.

#inter 0 0 lennard-jones $eps $sig $cut $shift 0


set tau_increase 10.

set delta_t [expr $integ_steps*$my_step]

puts "delta_t is, $delta_t"

set cap_ini 20.
set cap_fin 200.
set delta_cap [expr $cap_fin-$cap_ini]
set lin_coeff [expr $delta_cap/$tau_increase]

set n_iter_cap [expr $tau_increase/$delta_t]

set n_iter_cap [expr round($n_iter_cap)]

puts "the number of iterations while ramping the potential is, $n_iter_cap"


puts "delta_cap is, $delta_cap"

set cap $cap_ini

analyze set chains 0 1 $n_part



#for {set i 0} { $i <= $n_iter_cap } { incr i} {
#puts "t=[setmd time] E=[analyze energy total]"
#inter ljforcecap $cap; integrate $integ_steps

#set cap  [expr $cap + $lin_coeff*$delta_t ]

#set min [analyze mindist]

#}

#puts "Warmup finished. Minimal distance now $min"
#uncap forces
inter ljforcecap 0

puts "end of equilibration"

set vmd "no"

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

# The additional command imd steers the socket connection   #
# to VMD, e.g. sending the actual coordinates               #
    imd positions
}


#prepare the saving of the results




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"] 





#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 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 $obs3 "[setmd time] $rg"

puts $obs4 "[setmd time] [part 0 print pos]"

puts $obs5 "[setmd time]"





#puts "save configuration, [part ]"
part 0 print pos
#set f [open "config_$i" "w"]
#blockfile $f write tclvariable {box_l density}
#blockfile $f write variable box_l
#blockfile $f write particles {pos}
#close $f


set f [open "config_vel_$i" "w"]
#blockfile $f write tclvariable {box_l density}
#blockfile $f write variable box_l
blockfile $f write particles "id pos v"
close $f

integrate $integ_steps


# if you have turned on the vmd option you can now
    # follow what your simulation is doing
    if { $vmd == "yes" } { imd positions }


}

#close $obs


close $obs
close $obs2
close $obs3
close $obs4
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"