[ESPResSo] Re: Radius of Gyration

[Lorenzo Isella]

Hello,
You are right, I missed the setting up of a chain made up of all the 
particles I have.


> See the user's guide: analyze rg requires you to specify the topology, either 
> directly or before via "analyze set chains". In your case, you don't have a 
> real chain, but you can just tell it to treat your particles as one with e.g.
> set rg_temp [analyze  rg 0 1 $n_part]

However I am puzzled: I am running some test simulations with a small 
number of particles interacting via a very deep Lennard-Jones (LJ) 
potential. The simulations run within minutes on my box and I finally 
end up with a single aggregate constituted by the all the particles I 
have in the box. After the interaction has given rise to this aggregate, 
it is meaningful to calculate its radius of gyration (which I anyway 
start calculating ).
I tried first a simulation with 10 particles and then with 40 particles 
in a box of size 20 (first line of the code given in the following).
I must be either setting up poorly my system or misunderstanding 
something about how the radius of gyration is calculated by Espresso, 
since I have a larger value for a system with 10 rather than 40 
particles...and the interaction is the same and the box size is the same.
The code I re-paste here is quite close to one of the Espresso 
tutorials; one simply sets n_part to 10 or 40 and sees what the radius 
of gyration does.
Any hint at what is going wrong?
Many thanks

Lorenzo

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

set box_l 20.
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]]

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



puts "Ok after allocating a few particles"




set tot_time 1000.


set my_step 0.005


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 1000 

#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 2. ; set gamma  1.
thermostat langevin $temp $gamma





set sig 1.0
set cut 1.1
set eps 20.
set shift 0.
#set e1 12.
#set e2 10.
set r_off 0.
set omega 0.3

#inter 0 0 lj-cos2 $eps $sig $cut $shift $r_off $e1 $e2

inter 0 0 lj-cos2 $eps $sig $r_off $omega




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




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"


analyze set chains 0 1 $n_part

set rg_temp [analyze  rg 0 1 $n_part ]
puts "the radius of gyration is, $rg_temp"


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




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


}

incr count

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 }


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

}

#close $obs
close $obs
close $obs2
close $obs3

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"


--
Always remember to pillage before you burn.