#################################
####    PARAMETERS       ########
#################################
#temperature of the system
variable Temp equal 1.0

variable tdamp equal 1.0

#Equilibration
variable run1 equal 10000

#Restart every these timesteps
variable trestart equal 10000

#Dump data every these timesteps
variable dumpfreq equal 100

#Dump trajectory every these timesteps
variable dumptraj equal 100

variable seedthermo equal 125485

################################
####	DEFINTIONS 	    ########
################################
#Where the restart files of the run (non the equilibration) will be stored
variable folder index DUMPS
shell    mkdir ${folder}

#The starting configuration is always the same. Also the name of the restart files
variable rname   index Eq.5_1.N1000.data
variable simname index Run.R

units lj
atom_style angle 
boundary        p p p
neighbor 1.2 bin
neigh_modify every 1 delay 1 check yes 

restart ${trestart} ${folder}/Restart.${simname}.
read_data ${rname}


####################################
####	ANGLE	             #######
####################################
angle_style   cosine
angle_coeff  1 5.0 #20 sigma for realistic DNA (1 bead = 2.5 nm)

#######################################################
####	PAIRS -- REPULSIVE + TOPO II 	#######
#######################################################
pair_style  hybrid lj/cut 1.12246 soft 1.12246
pair_modify     shift yes
pair_coeff      * * lj/cut 1.0 1.0 1.12246
## NOW THERE IS TOPO2 -- TYPE 4 BEADS ARE SOFT AND CAN BE CROSSED
pair_coeff      * 4 soft 2.0  1.12246
pair_coeff      * 5 soft 20.0 1.12246

###################################
####  FENE	 ##################
##################################
bond_style   fene
#It turns off the LJ part of the FENE-bond potential between 1-2 atoms. However, the pairwise LJ interaction is still on
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
bond_coeff    2  30.0   1.6  1.0  1.0

####################################
####	THERMO AND INTEGRATION
####################################
thermo 1000
thermo_style   custom   step  temp  epair vol cpu
timestep 0.01

####################################
####	DUMPS	             #######
####################################
dump 2 all custom ${dumpfreq} ${folder}/${simname}.* id type mol mass x y z ix iy iz vx vy vz fx fy fz
dump_modify 2 sort id

##XYZ
dump 4 all xyz ${dumptraj} ${folder}/traj.xyz
dump_modify 4 element O N Ca P Au

compute 	MSDs all msd
compute		GYRs all gyration
fix imprimir all print ${dumpfreq} "$(step)  $(c_GYRs) $(c_MSDs[4])" screen no append radius_of_gyration.txt

#Compute total energies
compute ekin all ke
compute epot all pe
variable ekin equal c_ekin
variable epot equal c_epot
variable etot equal c_ekin+c_epot
fix energia all print ${dumpfreq} "$(step) ${ekin} ${epot} ${etot}" screen no append energy_total.txt

####################################
####    FIXES                #######
####################################
fix 1 all  nve
fix 2 all  langevin   ${Temp} ${Temp} ${tdamp}  ${seedthermo}

run ${run1}