Dear Kai,
I'm trying with ext_force and I see the chain is nicely stretched but still can't read it from numbers.
D.
######################################################################################
# #
# START TIME #
# #
######################################################################################
set systemTime0 [clock seconds]
######################################################################################
# #
# RANDOM SEED #
# #
######################################################################################
set n_cpu 8
set ran_seed [clock seconds]
set ran_seed [expr $ran_seed - 1407918000]
expr srand($ran_seed)
set cmd "t_random seed"
for {set i 0} {$i < $n_cpu} { incr i } { lappend cmd [expr $ran_seed+$i] }
eval $cmd
puts "$ran_seed"
######################################################################################
# #
# ESPRESSO FEATURES #
# #
######################################################################################
require_feature LENNARD_JONES
require_feature EXTERNAL_FORCES
require_feature PARTIAL_PERIODIC
require_feature LANGEVIN_PER_PARTICLE
require_feature MASS
require_feature EXCLUSION
require_feature BOND_CONSTRAINT
require_feature BOND_ANGLE_HARMONIC
require_feature OLD_DIHEDRAL ;# harmonic
#require_feature GAUSSIAN
######################################################################################
# #
# DEFAULTS #
# #
######################################################################################
set pbc_type "part_periodic" ;# periodic or non-periodic?
set constrained "yes" ;# do constrains exist or not?
set PI 3.1415926535
set accur 0.000001 ;# accuracy
######################################################################################
# #
# TOPOLOGY #
# #
######################################################################################
set NPART 100
set NBOND 99 ;# for circle NBOND=NPART
set vtffilename "measure.vtf" ;# trajectory file
set forcefilename "force.txt" ;# trajectory file
######################################################################################
# #
# SIMULATION #
# #
######################################################################################
set LX [expr 300]
set temperature 2.5 ;# simulation temperature
set gammat 1.0 ;# friction parameter Langevin thrmstt 1.0
set gammar 1.0 ;# friction parameter reel particles (3.0) 1.0
set massr 1.0 ;# mass reel particles
set timestep 0.05 ;# timestep
set num_frames 2000 ;# simulation cycles <================= THIS CHANGES =<<
set steps_per_frame 1000 ;# 1e4 : iterations per cycle <================= THIS CHANGES =<<
set frames_to_cnf 1 ;# 10: collect data each XXth cycle <================= THIS CHANGES =<<
set do_warmup "yes" ;# include warmup stage
set do_md "yes" ;# do MD
set save_warmingtrj "yes" ;# save conformational evolution trajectory, yes or no?
set warmup_maxenergy [expr 12.0*$NPART] ;# obsolete
set warmup_cap_min 10.0 ;# 10 : energy cap starts from this value
set warmup_cap_step 10.0 ;# 10 : liberate cap every XX steps
set warmup_steps 300 ;# 100 : thermalization integrations <================= THIS CHANGES =<<
set warmup_cycles 100 ;# 100 : thermalization cycles <================= THIS CHANGES =<<
set startsim 0 ;# 0: starting frame if not restarted
######################################################################################
# #
# SIMBOX #
# #
######################################################################################
#setmd box_l $LX $LX $LX
setmd box_l 300.0 300.0 300.0
setmd time_step $timestep
setmd skin [expr 15.]
setmd periodicity 1 1 1
######################################################################################
# #
# INTERACTIONS #
# #
######################################################################################
# Non-bonded interactions params
set ljsigma 1.0 ;# LJ sigma
set ljepsilon 1.0 ;# LJ epsilon
set ljrcut [expr $ljsigma * pow (2., 1./6.)] ;# LJ cut off
set ljcshift [expr 0.25 * $ljepsilon] ;# LJ shift
set ljroff 0.0
set ljrcap 3.0
set ljrmin 0.0
inter 0 0 lennard-jones $ljepsilon $ljsigma $ljrcut $ljcshift $ljroff $ljrcap $ljrmin ;# real-real
# Bonded interactions params HARMONIC
inter 1 harmonic 10. 1.0
# Bonded interactions params BENDING
set stiffness 100.0 ;# Stiffness
inter 2 angle $stiffness
######################################################################################
# #
# PARTICLES & INTERACTIONS #
# #
######################################################################################
#CREATING CHAIN
for {set pid 1} {$pid <= $NPART} {incr pid} {
set x [expr $pid]
set y [expr 0.0]
set z [expr 0.0]
part $pid pos $x $y $z mass $massr gamma $gammar
}
#FIXING ONE END OF CHAIN IN SPACE
part 1 fix 1 1 1
#INTRODUCING BONDED FENE
for {set pid 2} {$pid <= $NPART} {incr pid} {
part [expr $pid-1] bond 1 $pid
}
#INTRODUCING BENDING
set rid 1
set bid 2
for {set pid 3} {$pid <= $NPART} {incr pid} {
part [expr $pid-1+($rid-1)*$NPART] bond $bid [expr $pid-2+($rid-1)*$NPART] [expr $pid+($rid-1)*$NPART]
}
######################################################################################
# #
# OUTPUT FILES #
# #
######################################################################################
#WRITE CONFIGURATIONS
set vtffile [open $vtffilename w]
puts $vtffile "atom 0:[expr $NPART-1] radius 1 name DNA"
puts $vtffile "bond 0::[expr $NPART-1]"
puts $vtffile "unitcell $LX $LX $LX 90 90 90"
writevcf $vtffile
flush $vtffile
#WRITE FORCES
set forcefile [open $forcefilename w]
######################################################################################
# #
# THERMOSTAT #
# #
######################################################################################
thermostat langevin $gammat $temperature
puts "thermostat=[thermostat]"
puts "oki"
######################################################################################
# #
# WARM-UP #
# #
######################################################################################
if {$TPLG!="restart"} then {
puts "Warming up..."
set warmup_maxenergy [expr $NPART] ;# obsolete
set warmup_cap_min 1.0 ;# 10 : energy cap starts from this value
set warmup_cap_step 1.0 ;# 10 : liberate cap every XX steps
set warmup_steps 1000 ;# 100 : thermalization integrations <================= THIS CHANGES =<<
set warmup_cycles 5 ;# 100 : thermalization cycles <================= THIS CHANGES =<<
set save_warmingtrj "yes" ;# save conformational evolution trajectory, yes or no?
set wcap $warmup_cap_min
set i 1
if {$do_warmup=="yes"} then {
while {$i < $warmup_cycles} {
if {$i< 50} then {
setmd time_step 0.001
} else {
setmd time_step $timestep
}
set i [expr $i+1]
set wcap [expr $wcap + 10]
inter forcecap $wcap
integrate $warmup_steps
if { $save_warmingtrj == yes } then {
# writevcf $vtffile
# flush $vtffile
puts "[part 1 print]"
}
set e_tot [analyze energy total]
puts -nonewline [format "|Warmup | Step: %4d Energy: %.2f" $i $e_tot]\r
flush stdout
}
}
puts ""
puts "|Warmup | Finished. "
inter forcecap 0
}
######################################################################################
# #
# MAIN MD SIMULATION & PRODUCTION RUN #
# #
######################################################################################
#moto location
puts ""
puts "Main simulation started.."
if {$do_md=="yes"} then {
for {set frame 1} {$frame <= 10000000} {incr frame} {
integrate 1000
puts -nonewline [format "|Simulation | Step: %4d" $frame ]\r
flush stdout
writevcf $vtffile
flush $vtffile
part 100 ext_force [expr 50.0*sin($frame/2.0*$PI/180)] 0.0 0.0
set f(1) [lindex [part 98 print f] 0]
set f(2) [lindex [part 98 print f] 1]
set f(3) [lindex [part 98 print f] 2]
set Ft [expr pow($f(1)*$f(1)+$f(2)*$f(2)+$f(3)*$f(3),0.5)]
puts "Total force: $Ft Force vector: $f(1) $f(2) $f(3) Applied force: [expr 50.0*sin($frame/2.0/$PI/100)]"
puts $forcefile "$Ft $f(1) $f(2) $f(3) [expr 50.0*sin($frame/2.0/$PI/100)]"
}
}
close $forcefile
close $vtffile
######################################################################################
# #
# EXIT SIMULATION #
# #
######################################################################################
set systemTime1 [clock seconds]
set dtime [expr $systemTime1-$systemTime0]
puts "Elapsed in $dtime seconds."