[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 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
# 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.
- [ESPResSo] Re: Radius of Gyration,
Lorenzo Isella <=