#####################################################
#                                                   #
# Execute the script through:                       #
# lmp_foo < *.txt                                   #
#                                                   #
#####################################################



# VARIABLES
variable    T  equal 1.0
variable fx atom fx
variable fy atom fy
variable fz atom fz
variable vx atom vx
variable vy atom vy
variable vz atom vz
variable rx atom x
variable ry atom y
variable rz atom z

# Initialization
units		lj
boundary	p p p
molecule	scAB solvoChain_l20.txt
atom_style      template scAB
region 		myBox block 0 100 0 100 0 100
create_box	2 myBox #bond/types 2
create_atoms	0 random 216 98734 myBox mol scAB 4852

print "neighbor lists"
neighbor	0.7 bin
neigh_modify every 1 delay 0 check no
print "Set interactions"
velocity all create ${T}  9822 #4674 #7654
# pair potentials
pair_style morse 5.000000 #cut-off 10e-8
pair_coeff 1 1    0.71400    1.8000    1.0000    4.466
pair_coeff 1 2    0.01785    6.0000    1.7500    2.372
pair_coeff 2 2    0.01785    6.0000    2.5000    3.122

pair_modify shift yes
bond_style fene/expand
#CAUTION HERE bonds fene ignores morse interaction between connnected monomers
special_bonds lj 1.0 1.0 1.0 #gives exact energy with MC
##special_bonds fene
bond_coeff 1   120.0000    0.6375    0.0000    0.0000    3.1875
bond_coeff 2    60.0000    0.7125    0.0000    0.0000    3.5625


compute RGT all gyration/molecule
variable avRGT equal ave(c_RGT)
#compute bonds all bond/local dist
#variable maxBond equal max(c_RGT)
#compute 1 all property/local batom1 batom2 
#compute 2 all bond/local dist 
#compute prob all pe/atom bond
compute CM all com/molecule 
compute CMFx all atom/molecule v_fx
compute CMFy all atom/molecule v_fy
compute CMFz all atom/molecule v_fz
compute CMVx all atom/molecule v_vx
compute CMVy all atom/molecule v_vy
compute CMVz all atom/molecule v_vz
compute PE all pe/atom
print "Set simulation details"

fix therm all nvt temp $T $T 0.09
fix_modify therm energy yes
#fix aveRGT all ave/time 1000 1 1000 c_RGT mode vector file out_RGT.dat
#fix aveCM all ave/time 1000 1 1000 c_CM[1] c_CM[2] c_CM[3] c_CMFx c_CMFy c_CMFz file output/com.* mode vector title1 "cm_x cm_y cm_z cm_fx cm_fy cm_fz"
#fix volumeChange all deform 1 x final 0.0 20.0 y final 0.0 20.0  z final 0.0 20.0 units box
##############################################################
thermo_style custom step temp pe ke epair ebond etotal vol fmax v_avRGT #v_maxBond 
thermo_modify norm no flush yes lost error
thermo 2000
restart 2500000 conf_restart.*


run	0
variable t equal pe
print "Initial potential energy = $t"

dump dumpWrap all xyz 200 out_dump.xyz
#dump dumpForces all custom 1 forces.dat id fx fy fz c_PE
#dump forces all custom 1 forces.dat v_rx v_ry v_rz v_fx v_fy v_fz
#dump 1 all local 10 tmp.dump index c_1[1] c_1[2] c_2
timestep 0.0005

run 250000000


print "All done"