#NP_A_Diamondene


variable       T2 equal 350
variable       T1 equal 340
variable       y1 equal 0
variable       y2 equal 0
timestep       0.001




##########################################
##########################################
label           loop
variable        i loop 2 2000

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

units          metal
dimension 	    3 
boundary       p p s
atom_style 	molecular


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

if "$i ==1 " then &
    "read_data 	NP_A_Diamondene.data" &
    "variable       stemp equal 1" &
elif "$i ==2 " &
"variable       stemp equal 345" &
"read_data 	relaxed_NP_A_Diamondene.data" &
else "read_data 	relaxed_NP_A_Diamondene.data"

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

group G1 type 1
group G2 type 2
group G3 type 3
group G4 type 4


##---------------FORCE FIELDS---------------------------------
bond_style 	harmonic
bond_coeff * 0.0 1.0
special_bonds lj/coul 1.0 1.0 1.0


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

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

neighbor       2.0 bin
neigh_modify   every 1


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

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


##---------------SETTINGS-------------------------------------
reset_timestep 0
compute            pePerAtom all pe/atom 
compute            Peall all pe 
compute            2 all property/local batom1 batom2 btype
compute            3 all bond/local dist


velocity       all create ${stemp} 12345 units box dist gaussian 
fix            1 all npt temp ${stemp} ${stemp} 0.1 x 0 0 20 y 0 0 20 drag 0.1 
fix            2 all momentum 1 linear 1 1 1 angular rescale


compute            11234 G1 coord/atom cutoff 2.0 1 2 3 4
compute            2234 G2 coord/atom cutoff 2.0 2 3 4
compute            334 G3 coord/atom cutoff 2.0 3 4
compute            44 G4 coord/atom cutoff 2.0 4 

compute            BondsSum all reduce sum c_11234[2] c_2234[2] c_334[2] c_11234[1] &
c_2234[1] c_334[1] c_44 c_11234[4] c_11234[3] c_2234[3]

variable            inBrokenBonds equal 324-c_BondsSum[1]
variable            midBrokenBonds equal 108-c_BondsSum[2]
variable            outBrokenBonds equal 324-c_BondsSum[3]

variable            11Bonds equal c_BondsSum[4]*0.5
variable            22Bonds equal c_BondsSum[5]*0.5
variable            33Bonds equal c_BondsSum[6]*0.5
variable            44Bonds equal c_BondsSum[7]*0.5

variable            totalBrokenBonds equal v_inBrokenBonds+v_midBrokenBonds+v_outBrokenBonds
variable            recreatedBonds equal &
(c_BondsSum[4]+c_BondsSum[5]+c_BondsSum[6]+c_BondsSum[7])*0.5+c_BondsSum[8]+c_BondsSum[9]+c_BondsSum[10]
variable            grossBrokenBonds equal v_totalBrokenBonds-v_recreatedBonds

if " $(abs(v_T2-v_T1)) <70" then &
	"fix            9 all halt 400 v_totalBrokenBonds > 1.000000 error continue"

#In order to avoid errors czed from blown away atoms, deleted and leaving bonds behind.
# refresh/timestep bondtype Rmax probablity(100% seed)
fix 3 all bond/break 1 1 10.0 prob 1 85784
fix 4 all bond/break 1 2 10.0 prob 1 85784
fix 5 all bond/break 1 3 10.0 prob 1 85784
#fix 6 all bond/break 1 4 10.0 prob 1 85784
#fix 7 all bond/break 1 5 10.0 prob 1 85784

##---------------OUTPUT--------------------------------------
thermo_style   custom step temp v_grossBrokenBonds v_totalBrokenBonds v_inBrokenBonds v_midBrokenBonds &
v_outBrokenBonds v_11Bonds v_22Bonds v_33Bonds v_44Bonds c_BondsSum[8] c_BondsSum[9] c_BondsSum[10]
dump           ${i} all custom 1000 NP_A_DNT_6_${stemp}C id type x y z c_pePerAtom
thermo            200
thermo_modify   lost warn
thermo_modify  format 6 %4.0f
thermo_modify   format 7 %4.0f
thermo_modify   format 8 %6.0f
thermo_modify  format 9  %3.0f
thermo_modify   format 10 %3.0f
thermo_modify   format 11 %3.0f
thermo_modify   format 12 %5.0f
thermo_modify   format 13 %3.0f
thermo_modify   format 14 %3.0f

##---------------RELAXATION--------------------------------------
if "$i ==1 " then &
"run 500" &
else    "run 500" 

##---------------BiSection method--------------------------------------
if " $(abs(v_T2-v_T1)) >70" then &
	"""if "${grossBrokenBonds}<1.000000 && $i>1" then 
		'if "${T1}>${T2}" then 
			"variable       T1 equal ${stemp}" 
			"variable       stemp equal ${T1}+400" 
		else "variable       T1 equal ${stemp}" 
		     "variable       stemp equal (v_T1+v_T2)/2" ' 
	else	"variable       T2 equal ${stemp}" 
		"variable       stemp equal (v_T1+v_T2)/2" """ & 
elif " $(abs(v_T2-v_T1)) >1" &
	"""if "${grossBrokenBonds}<1.000000" then 
		"variable       y1 equal v_grossBrokenBonds-1.000000"
		"variable       T1 equal ${stemp}" 
		"variable       stemp equal v_T1+v_y1 * (v_T2-v_T1)/(v_y2-v_y1)" 
	else	"variable       T2 equal ${stemp}" 
		"variable       y2 equal v_grossBrokenBonds-1.000000"
		"variable       stemp equal v_T2-v_y2 * (v_T2-v_T1)/(v_y2-v_y1)" """ &
else "quit" 


if "$i==1" then &
    "write_data 	relaxed_NP_A_Diamondene.data nocoeff nofix"


clear
next            i
jump            NP_A_Diamondene.in loop

##########################################
##########################################