[ESPResSo] Question about bond-angle interaction

[Lorenzo Isella]

Dear All,
I have already posted on this list, mainly to ask whether I could used
Espresso for simulating the dynamics of nanoparticles rather than
molecules/polimers.
My problem is to have these particles "sticking" when they collide. This
is essential for me, since I want to study the complicated, absolutely
non-spherical, aggregates they give rise to.
I have been toying around both with Lennard-Jones and Morse potential.
The problem (which I think is due to the potentials being a function of
the distance between two particles only) is that an aggregate is indeed
formed, but I always end up with a compact, almost spherical structure.
To make it less an abstract discussion, here is the code I have used:

set n_part 20; set density 0.07
set box_l [expr pow($n_part/$density,1./3.)]

setmd box_l $box_l $box_l $box_l
setmd periodic 1 1 1



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
}

setmd time_step 0.01; setmd skin 0.4
set temp 0.0005; set gamma 0.01
thermostat langevin $temp $gamma





#set sig 1.0; set cut [expr 1.12246*$sig]
#set eps 1.0; set shift [expr 0.25*$eps]
#inter 0 0 lennard-jones $eps $sig $cut $shift 0
inter 0 0  morse 1000. 0.01 1. 1.3

#inter 3 angle 1000.0 3.14



set integ_steps 200

for {set cap 20} {$cap < 200} {incr cap 20} {
puts "t=[setmd time] E=[analyze energy total]"
inter ljforcecap $cap; integrate $integ_steps
set min [analyze mindist]

}

puts "Warmup finished. Minimal distance now $min"
#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 obs [open "rg_lorenzo.dat" "w"]

set integ_step 20

for {set i 0} { $i < 1000 } { incr i} {
set temp [expr [analyze energy kinetic]/(1.5*$n_part)]
puts "t=[setmd time] E=[analyze energy total], T=$temp"
integrate $integ_steps

puts "save configuration"

set f [open "config_$i" "w"]
blockfile $f write tclvariable {box_l density}
blockfile $f write variable box_l
blockfile $f write particles {id pos type}
close $f
# if you have turned on the vmd option you can now
   # follow what your simulation is doing
   if { $vmd == "yes" } { imd positions }

#Now I calculate the radius of gyration
# set rg [lindex [analyze rg] 0]
#puts $obs "[setmd time] $rg"

}

#close $obs
puts "end of integration"

#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]
}


Any suggestions about a different interaction to use? In a previous post
I was suggested a combination of Lennard-Jones and a stochastically
added FENE bond, but this is still obscure to me.
So, I was thinking about some interaction which depends on the angle
determined by 3 particles (bond-angle interaction). However, setting the
interaction in the previous code as:

inter 4 angle 1000.0 3.14
part 0  bond 4 0 0

Leads to some strange errors. Actually, the number after part should
probably be the individual number stuck on each particles and not the
type but then how can I say that I want all the particles to interact
according to that potential? And then:   "angle" stands for
BOND_ANGLE_COSINE  according to my myconfig.h file.

I am a bit confused here, I may be using funny numbers in the code. Any
help is really appreciated.
Many thanks

Lorenzo