#	0 --- 0 Harm
#	|     |		+ DIHEADRAL
#	0 --- 0	Harm
#	|     |
#	0 --- 0

#variable rname index trefoil_lammps.out 
variable rname index Trefoil.Patchy.Restart.20000
variable simname index Trefoil.Patchy

units lj

atom_style molecular #with stiffness and dihedral (torsion)

boundary        p p p

###
# Pair interactions require lists of neighbours to be calculated
###
#neighbor 0.3 bin
neighbor 1.9 bin
neigh_modify every 1 delay 1 check yes 

### Generate RESTART file, SPECIAL format, not a .txt file ###
### Useful if simulation crashes ###
restart 10000 ${simname}.Restart
#############################

#processors 2 1 1

### READ "start" data file ###
#read_data ${rname}
######################

### READ "REstart" data file ###
read_restart ${rname}
######################

### reset timestep ###
#reset_timestep 0 
######################

###
# Define groups 
###
group backbone type 1
group patch type 2 
group all type *

######COMPUTE STUFF#######################################################################################3
compute 	Stress all stress/atom dihedral  angle
compute 	Curv all stress/atom angle
compute 	Twist all stress/atom dihedral
compute 	Bond all stress/atom bond

variable	BondTot atom sqrt((c_Bond[1]+c_Bond[2]+c_Bond[3]+c_Bond[4]+c_Bond[5]+c_Bond[6])^2)
variable	StressTot atom sqrt((c_Stress[1]+c_Stress[2]+c_Stress[3]+c_Stress[4]+c_Stress[5]+c_Stress[6])^2)

variable	CurvTot atom sqrt((c_Curv[1]+c_Curv[2]+c_Curv[3]+c_Curv[4]+c_Curv[5]+c_Curv[6])^2)
variable	TwistTot atom sqrt((c_Twist[1]+c_Twist[2]+c_Twist[3]+c_Twist[4]+c_Twist[5]+c_Twist[6])^2)

###DUMPS ##################################################################################################
# Dump configs + box coords
dump   1   all   cfg   100   movie/${simname}.0*.cfg  id  type xs ys zs id type v_StressTot v_CurvTot v_TwistTot v_BondTot
dump_modify 1 element O C
#dump   2 all atom 1000 dump4xmovie1
######################################################################################3
angle_style   hybrid harmonic cosine
angle_coeff   1 cosine 10.0
angle_coeff   2 harmonic 200.0 90.0

dihedral_style charmm	#K(1+cos(n phi - d))
dihedral_coeff 1 100.0 1 120 0.0 

#PAIR INTERACTION LJ#########################################################################
pair_style	lj/cut 1.12246152962189
pair_coeff	1 1 1.0 1.0 1.12246152962189
pair_coeff	1 2 0.0 1.0 1.12246152962189
pair_coeff	2 2 0.0 1.0 1.12246152962189

####################################################################################################
# Bond interaction between bonded atoms
#
# Fene potential + Lennard Jones 12-6:
#
#  E= - 0.5 K R0^2 ln[ 1- (r/R0)^2] 
#     + 4epsilon[ (sigma/r)^12 - (sigma/r)^6] + epsilon
#
bond_style   fene
#
# For style fene, specify:
#   * bond type
#   * K (energy/distance^2)
#   * R0 (distance)
#   * epsilon (energy)  (LJ component)
#   * sigma (distance)  (LJ component)
special_bonds fene#<=== I M P O R T A N T (new command)
bond_coeff   1  30.0   1.6  1.0   1.0
#######HYBRID BONDS##########################################################################################
#bond_style hybrid harmonic fene
#bond_coeff 1 harmonic 250.0 1.1
#bond_coeff 2 harmonic 250.0 1.1
#bond_coeff 2 fene 5.0 1.7 1.0 1.0 
#########################################################################################################
###
# Set up fixes
###

fix 1 all rigid molecule   
fix 2 all   langevin   1.0 1.0  2.0   424242   ###Langevin integrator Tstart Tstop 1/friction rndseed

compute 	MSDs backbone msd
#compute 	COM pulled com
compute		GYRs backbone gyration
#compute		GYR pulled gyration/molecule 
#compute		GYRT pulled gyration/molecule tensor

###STRESS########################################################################
compute		p all reduce sum c_Stress[1] c_Stress[2] c_Stress[3]
variable	press equal -(c_p[1]+c_p[2]+c_p[3])/(3*vol)
#####################################################################################
##### Sample thermodynamic info  (temperature, energy, pressure, etc.) #####
thermo 1000
thermo_style   custom   step  temp  pe  c_MSDs[4] c_GYRs press
###
# set timestep of integrator
###
timestep 0.01 #0.012

###
# run integration for a number of steps
###
run 20000