#DP_A_DNT_16_900C

##---------------INITIALIZATION-------------------------------

units          metal
dimension 	    3 
boundary       p s s
atom_style 	molecular

##---------------ATOM DEFINITION------------------------------

read_data 	DP_A_DNT_16.data


##---------------FORCE FIELDS---------------------------------
bond_style harmonic
bond_coeff * 0.0 1.0

pair_style 	airebo 3.0
pair_coeff     * * CH.airebo H C C C C H


special_bonds lj/coul 1.0 1.0 1.0



##---------------neighbor list---------------------------------

neighbor       2.0 bin
neigh_modify   every 1


##---------------minimize---------------------------------

min_style      sd
minimize       1.0e-4 1.0e-6 100 1000


##---------------RELAXATION--------------------------------------

variable       T2 equal 1
variable       T1 equal 900
timestep       0.001

variable       stemp equal 900
label           loop
variable        i loop 2000


velocity       all create ${stemp} 12345 units box dist gaussian 
fix            1 all npt temp ${stemp} ${stemp} 0.1 x 0 0 0.1 drag 0.1 
fix            2 all momentum 1 linear 1 1 1 angular rescale
# refresh/timestep bondtype Rmax probablity(100% seed)
fix 3 all bond/break 1 1 2.0 prob 1 85784
fix 4 all bond/break 1 2 2.0 prob 1 85784
fix 5 all bond/break 1 3 2.0 prob 1 85784
fix 6 all bond/break 1 4 2.0 prob 1 85784
fix 7 all bond/break 1 5 2.0 prob 1 85784


variable bondsNumber equal bonds
variable inBrokenBonds equal f_4[2]
variable interBrokenBonds equal f_5[2]
variable outBrokenBonds equal f_6[2]
variable totalBrokenBonds equal ${inBrokenBonds}+${interBrokenBonds}+${outBrokenBonds}
compute            1 all pe/atom 
compute            Peall all pe 

run             30

##---------------OUTPUT--------------------------------------

dump           1 all custom 1000 DP_A_DNT_16_900C id type x y z c_1

thermo         10
thermo_modify   lost warn
thermo_style custom step temp pe ke etotal v_inBrokenBonds &
v_interBrokenBonds v_outBrokenBonds v_totalBrokenBonds 


##---------------BiSection method-------------------------------------
if "${totalBrokenBonds}==0" && "${T1}>${T2}" then &
	"print 'the system is still stable because the broken bonds =${totalBrokenBonds}'" &
	"variable       T1 equal ${stemp}" &
	"variable       stemp equal ${T1}+400 " &
elif "${totalBrokenBonds}==0" && "${T1}<${T2}" &
	"variable       T1 equal ${stemp}" &
	"variable       stemp equal (${T1}+${T2})/2" &
elif "${totalBrokenBonds}>6" &
	"variable       T2 equal ${stemp}" &
	"variable       stemp equal (${T1}+${T2})/2" &
elif "${totalBrokenBonds}>0" && ${totalBrokenBonds}<6"  &
	"print 'The critical Temperature is = ${stemp}'" &
	"jump            2DDTube.in break" 

next            i
jump            2DDTube.in loop

label           break


#run            50000


#if you insist on dealing with this from inside of LAMMPS,
#you could define bonds and then use bond_style none
#(or harmonic with K=0 if that doesn't work). you only
#have to make sure you set you special_bonds settings
#to 1.0 1.0 1.0 so that no interactions are excluded.
#special_bonds lj/coul 1.0 1.0 1.0