espressomd-users
[Top][All Lists]

## [ESPResSo] Re: Radius of Gyration

 From: Lorenzo Isella Subject: [ESPResSo] Re: Radius of Gyration Date: Sat, 08 Dec 2007 21:35:30 +0100 User-agent: Icedove 1.5.0.14pre (X11/20071020)

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

```