Hello Georg Rempfer
This is my simulation script.
set n_solute 100
set n_solvent 2500
set box_l 200
setmd box_l $box_l $box_l $box_l
setmd periodic 1 1 1
for {set i 0} { $i < 100 } {incr i} {
set posx [expr $box_l*[t_random]]
set posy [expr $box_l*[t_random]]
set posz [expr $box_l*[t_random]]
set vx [gauss_random]
set vy [gauss_random]
set vz [gauss_random]
part $i pos $posx $posy $posz type 0 v $vx $vy $vz q -25 mass 1
}
for {set i 100} { $i < 2600 } {incr i} {
set posx [expr $box_l*[t_random]]
set posy [expr $box_l*[t_random]]
set posz [expr $box_l*[t_random]]
set vx [gauss_random]
set vy [gauss_random]
set vz [gauss_random]
part $i pos $posx $posy $posz type 1 v $vx $vy $vz q 1 mass 1
}
setmd time_step 0.006
# 1 TIMESTEP = 6E-14
setmd skin 0.4
set temp 2.493
# 300K
set gamma 0.4
thermostat langevin $temp $gamma
# WCA POTENTIAL
#radius of solute=16nm and radius of solvent=1.0nm
set sig 1.0
set eps 1.0
set cut [expr 1.22426*$sig]
set shift [expr 0.25*$eps]
inter 0 0 lennard-jones $eps $sig $cut $shift 0
inter 1 0 lennard-jones $eps $sig $cut $shift 0
inter 1 1 lennard-jones $eps $sig $cut $shift 0
inter coulomb 0.718 p3m tunev2 accuracy 1e-3
# p3m VARIABLES LIKE NO. OF MESH POINTS,CHARGE ASSIGNMENT ORDER ETC. ARE CALCULATED BY THE FUNCTION TUNEV2
set inter_steps 100000
for {set cap 30} {$cap < 200} {incr cap 10} {
set temp [expr [analyze energy kinetic]/(1.5*[setmd n_part ])]
puts "t=[setmd time] E=[analyze energy total] T=$temp"
inter ljforcecap $cap;
integrate $inter_steps
}
inter ljforcecap 0
set g [open "analysis.data" "w"]
set n_part [expr ($n_solute + $n_solvent)]
for {set i 0} { $i < 3350 } {incr i} {
puts "step $i ftime=[setmd time] energy=[analyze energy total]"
puts "temp = [expr [analyze energy kinetic]/(1.5*[setmd n_part])]"
integrate 300
set f [open "config_$i" "w"]
blockfile $f write tclvariable {box_l}
blockfile $f write variable box_l
blockfile $f write particles {id pos type v f q mass}
set temp [expr [analyze energy kinetic]/(1.5*[setmd n_part ])]
puts $f " \ { energy [analyze energy total] $temp}"
puts $g " [analyze energy kinetic] [analyze pressure total] [analyze energy total] $temp"
close $f
}
close $g
And my equilibration steps goes like this
t=0.0 E=5020.1630531796645 T=1.018479456496
t=599.999999998393 E=18220.59709794523 T=5.810625670095469
t=1199.9999999989084 E=15816.115978703594 T=5.256242903696039
t=1800.0000000074576 E=34210.471905675084 T=9.991591589316466
t=2400.0000000026675 E=36777.37784842733 T=10.675911375103668
t=2999.9999999884794 E=42742.20641527184 T=11.999132235697983
t=3599.9999999742913 E=26673.72159004423 T=7.974431903864867
t=4199.999999967986 E=4779348.10622052 T=1226.6668964027733
t=4799.999999999272 E=57093.543095210924 T=15.791968482255356
Thanks
Uday