[ESPResSo] Radius of Gyration

[Lorenzo Isella]

Dear All,
I have a problem when dealing with a test case (small system with a
few particles). For some reason, the calculation of the radius of
gyration returns NaN at any time. This has nothing to do with my
previous mail (I have not installed any CVS version of Espresso).
Consider for instance the following code:

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.01


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"

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

I am just allocating a few particles in a small box letting them
interact via a Lennard-Jones potential.
Why is the calculation of the radius of gyration failing?
Kind Regards

Lorenzo