[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[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 (X11/20071020)

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


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

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

reply via email to

[Prev in Thread] Current Thread [Next in Thread]