Dear all,
I am trying to simulate a system with 2 and 3 particle molecules, however on giving the command [analyze set] all the molecule IDs are shown to be 0. How can I give different IDs to the molecules.
Also, is there any way to observe the aggregation of a system such that I get the distribution of the aggregate size.
My script is given below:
set particle1 100
set particle2 100
set particles [expr $particle1+$particle2]
set density 0.001
set volume [expr ($particle1+$particle2)/$density]
set box_length [expr pow($volume, 1.0/3.0)]
set temp 2.500
set t_step 0.01
set gamma 1
setmd box_l $box_length $box_length $box_length
setmd time_step $t_step
setmd skin 0.4
thermostat langevin $temp $gamma
inter 0 fene 10.0 2.0 1.0
inter 1 fene 9.0 1.9 1.0
set molid 0
for {set i 0} {$i<$particle1} {incr i 2} {
set x [expr rand()*$box_length]
set y [expr rand()*$box_length]
set z [expr rand()*$box_length]
part $i pos $x $y $z type 0 v 0 0 0 mass 15
part [expr $i+1] pos [expr $x+1] $y $z type 0 v 0 0 0 mass 15
part $i bond 0 [expr $i+1] mol $molid
incr molid
}
for {set j $i} {$j<$particles} {incr j 3} {
set x [expr rand()*$box_length]
set y [expr rand()*$box_length]
set z [expr rand()*$box_length]
part $j pos $x $y $z type 0 v 0 0 0 mass 15
part [expr $j+1] pos $x [expr $y+1] $z type 0 v 0 0 0 mass 14
part [expr $j+2] pos $x $y [expr $z+1] type 0 v 0 0 0 mass 15
part [expr $j] bond 1 [expr $j+1] mol [expr $molid+1]
part [expr $j+1] bond 1 [expr $j+2] mol [expr $molid+1]
incr molid
}
analyze set chains 0 [expr $particle1/2] 2
analyze set chains $particle1 [expr $particles/3] 3
analyze set topo_part_sync
puts "[analyze set]"
inter 0 0 lennard-jones 1 1 2.5 2.0 0
inter 1 1 lennard-jones 1 1 2.5 2.0 0
inter 1 0 lennard-jones 1 1 2.5 2.0 0
#Hence we get the following reduced units for epsilon*=1;sigma*=1; m*= 1---- t*=1 ps; T*=120.27K; den*=1.6605kg/(m^3); numden*=1 nm^-3; P*=16.605bar
set vmd "ye"
if {$vmd == "yes"} {
prepare_vmd_connection Simple_System_vmd 100
imd positions
}
set dtem 100
#set inter_steps 1000
for {set cap 0} {[analyze mindist] < 0.9 || $dtem > 0.5} {incr cap} {
inter forcecap $cap
integrate 5
set tem1 [expr [analyze energy kinetic]/(1.5*[setmd n_part])]
set dtem [expr abs($tem1 - 2.5)]
puts "$cap [setmd temperature] $dtem"
}
inter forcecap 0
set rd [open "rdf.dat" "w"]
set fi [open "analysis.dat" "w"]
set ag [open "aggregation.dat" "w"]
for {set k 0} {$k<1000} {incr k} {
integrate 10
puts $fi "Step $k ----t : [setmd time] T : [expr [analyze energy kinetic]/(1.5*[setmd n_part])] E : [analyze energy total] K.E. : [analyze energy kinetic] P : [analyze pressure total]"
if {$vmd == "yes"} {imd positions}
puts $ag "[analyze aggregation 1.2 0 0]"
}
puts $rd "[analyze rdf 0 0]"
close $rd
close $ag
Regards
Udit Dewan