myprojects-hg-reborn

Fork

set n_part 200; set density 0.7
set box_l [expr pow($n_part/$density,1./3.)]

setmd box_l $box_l $box_l $box_l
setmd periodic 1 1 1

set q 1; set type 0
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 q [expr -$q]; set type [expr 1-$type]
part $i pos $posx $posy $posz q $q type $type
}

setmd time_step 0.01; setmd skin 0.4
set temp 1; set gamma 1
thermostat langevin $temp $gamma





set sig 1.0; set cut [expr 1.12246*$sig]
set eps 1.0; set shift [expr 0.25*$eps]
inter 0 0 lennard-jones $eps $sig $cut $shift 0
inter 1 0 lennard-jones $eps $sig $cut $shift 0
inter 1 1 lennard-jones $eps $sig $cut $shift 0
inter coulomb 10.0 p3m tunev2 accuracy 1e-3 mesh 32


set integ_steps 200

for {set cap 20} {$cap < 200} {incr cap 20} {
puts "t=[setmd time] E=[analyze energy total]"
inter ljforcecap $cap; integrate $integ_steps
set min [analyze mindist]

}

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

puts "end of equilibration"

set vmd "yes"

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 "rg_lorenzo.dat" "w"] 


for {set i 0} { $i < 3 } { incr i} {
set temp [expr [analyze energy kinetic]/(1.5*$n_part)]
puts "t=[setmd time] E=[analyze energy total], T=$temp"
integrate $integ_steps

puts "save configuration"

set f [open "config_$i" "w"]
blockfile $f write tclvariable {box_l density}
blockfile $f write variable box_l
blockfile $f write particles {id pos type}
close $f
# if you have turned on the vmd option you can now
    # follow what your simulation is doing
    if { $vmd == "yes" } { imd positions }

#Now I calculate the radius of gyration
# set rg [lindex [analyze rg] 0]
#puts $obs "[setmd time] $rg"

}

#close $obs
puts "end of integration" 

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