#	PATCHY WITH FIX RIGID
#	SIDEWAYS BEADS
#	1 -- 1
#	2 -- 2 
#	 ....
variable rname index trefoil_lammps_patchy.out
#variable rname index Trefoil.StiffTorsion.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 (atom type 1 is group 'all')
###
group backbone type 1
group patch type 2 
group all type *

######STRESS######################3
compute 	Stress all stress/atom dihedral  angle
compute 	Curv all stress/atom angle
compute 	Twist all stress/atom dihedral

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)


#############################################
# 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
dump_modify 1 element O C
#dump   2 all atom 1000 dump4xmovie1

angle_style   hybrid harmonic cosine
angle_coeff   1 cosine 20.0
angle_coeff   2 harmonic 200.0 90.0

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

#SOFT/////////////////////////////////////////////
pair_style     soft 1.12246152962189
#
#  pair_coeff for soft, specify 3:
#    * atom type interacting with
#    * atom type 
#    * A
#    * cutoff
 pair_coeff 1 1  100.0  1.0
 pair_coeff 1 2  0.0  0.673477
 pair_coeff 2 2  0.0   0.224492

#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   3  1.0   1.0

#variable prefactor equal ramp(10,1.5)
#fix 4 rings adapt 1 bond fene r0 * * v_prefactor 
###################################################

###################################################################
# Pair interaction between bonded atoms
#
# Harmonic potential: 
#
#  E= K * (R-R0)^2
#

bond_style   harmonic 

#
# For style harmonic, specify:
#   * bond type
#   * K (energy/distance^2) 
#   * R0 (distance)
#

bond_coeff   1  100.0   1.1
#bond_coeff   2  250.0   1.1

#bond_coeff   2   harmonic 0.1   1.0 
###################################################

###
# Set up fixes
###

#fix 1 all nve
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
 
 
#SECOND RUN, UNFIX SOFT, USE LJ/////////////////////////


#unfix 4
#pair_style	lj/cut 1.12246152962189
#pair_coeff	* * 1.0 1.0 1.12246152962189
#pair_modify shift yes
#thermo		1000


#run		19000