#########################################################################################
#
# Diffusion of a single polymer with hydrodynamic interactions (Lattice-Boltzmann)
#
#########################################################################################
# length of the polymer chain
set chain_length 16

setmd box_l 32. 32. 32.
setmd skin 0.2
setmd time_step 0.01
cellsystem domain_decomposition -no_verlet_list

inter 0 0 lennard-jones 1. 1. 1.226 0.25 0.
inter 0 fene 7. 2.

set vtffile [ open "polymer-diffusion.vtf" "w" ]
set posfile [ open "pos.dat" "w" ]
set vfile [ open "vel.dat" "w" ]

polymer 1 $chain_length 1. bond 0 mode PSAW

writevsf $vtffile

puts "Warming up the polymer chain."
## For longer chains (>100) an extensive 
## warmup is neccessary ...
setmd time_step 0.002
thermostat langevin 1. 10.
for { set i 0 } { $i < 100 } { incr i } {
  inter ljforcecap $i
  integrate 1000
}
puts "Done."
inter ljforcecap 0.
integrate 10000
setmd time_step 0.01
integrate 50000

# remove the langevin thermostat
thermostat off
lbfluid grid 1. dens 1. visc 5. tau 0.1 friction 5.
thermostat lb 1.

puts "Warming up the system with LB fluid."
integrate 1000
puts "Done."

for { set i 0 } { $i < 1000 } { incr i } {
    integrate 100
    set vx 0.
    set vy 0.
    set vz 0.
    for { set j 0 } { $j < [ setmd n_part ] } { incr j } {
	set v [ part $j print v]
        set vx [ expr $vx + [ lindex $v 0 ] ]
        set vy [ expr $vy + [ lindex $v 1 ] ]
        set vz [ expr $vz + [ lindex $v 2 ] ]
    }
    writevcf $vtffile
    puts $posfile "[ setmd time ] [ analyze centermass 0 ]" 
    puts $vfile "[ setmd time ]  [ expr $vx/$chain_length ] [ expr $vy/ $chain_length ] [ expr $vz/ $chain_length ]"
}

close $vtffile
close $posfile
close $vfile