#DP_A_DNT_20 variable TotalTimeSteps equal 500000 variable T2 equal 1100 variable T1 equal 1650 timestep 0.001 ########################################## ########################################## label loop variable i loop 2 2000 ##---------------INITIALIZATION------------------------------- units metal dimension 3 boundary p s s atom_style molecular ##---------------RELAXATION-------------------------------------- if "$i ==1 " then & "read_data DP_A_DNT_20.data" & "variable stemp equal 1" & elif "$i ==2 " & "variable stemp equal 2050" & "read_data relaxed_DP_A_DNT_20.data" & else "read_data relaxed_DP_A_DNT_20.data" ##---------------ATOM DEFINITION------------------------------ group H_in type 1 group G1 type 2 group G2 type 3 group G3 type 4 group G4 type 5 group H_out type 6 ##---------------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 H C C C C H ##---------------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 velocity all create ${stemp} 12345 units box dist gaussian fix 1 all npt temp ${stemp} ${stemp} 0.1 x 0 0 20 drag 0.5 fix 2 all momentum 1 linear 1 1 1 angular rescale compute H1 H_in coord/atom cutoff 2.0 2 compute 11 G1 coord/atom cutoff 2.0 2 compute 12 G1 coord/atom cutoff 2.0 3 compute 13 G1 coord/atom cutoff 2.0 4 compute 14 G1 coord/atom cutoff 2.0 5 compute 22 G2 coord/atom cutoff 2.0 3 compute 23 G2 coord/atom cutoff 2.0 4 compute 24 G2 coord/atom cutoff 2.0 5 compute 33 G3 coord/atom cutoff 2.0 4 compute 34 G3 coord/atom cutoff 2.0 5 compute 44 G4 coord/atom cutoff 2.0 5 compute H2 H_out coord/atom cutoff 2.0 5 compute H_inBonds all reduce sum c_H1 compute H_outBonds all reduce sum c_H2 compute inBonds all reduce sum c_12 compute midBonds all reduce sum c_23 compute outBonds all reduce sum c_34 compute 11Bonds all reduce sum c_11 compute 22Bonds all reduce sum c_22 compute 33Bonds all reduce sum c_33 compute 44Bonds all reduce sum c_44 compute 14Bonds all reduce sum c_14 compute 13Bonds all reduce sum c_13 compute 24Bonds all reduce sum c_24 variable H_inBrokenBonds equal 360-c_H_inBonds variable H_outBrokenBonds equal 360-c_H_outBonds variable inBrokenBonds equal 1080-c_inBonds variable midBrokenBonds equal 360-c_midBonds variable outBrokenBonds equal 1080-c_outBonds variable 11Bonds equal c_11Bonds*0.5 variable 22Bonds equal c_22Bonds*0.5 variable 33Bonds equal c_33Bonds*0.5 variable 44Bonds equal c_44Bonds*0.5 variable totalBrokenBonds equal v_inBrokenBonds+v_midBrokenBonds+v_outBrokenBonds variable recreatedBonds equal (c_11Bonds+c_22Bonds+c_33Bonds+c_44Bonds)*0.5+c_14Bonds+c_13Bonds+c_24Bonds variable grossBrokenBonds equal v_totalBrokenBonds-v_recreatedBonds fix 9 all halt 400 v_totalBrokenBonds > 3.333333 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-------------------------------------- variable 1 equal 0.2/(tpcpu+0.0000000001) thermo_style custom step temp pe ke & v_grossBrokenBonds v_totalBrokenBonds v_inBrokenBonds v_midBrokenBonds & v_outBrokenBonds v_11Bonds v_22Bonds v_33Bonds v_44Bonds c_14Bonds c_13Bonds c_24Bonds v_H_inBrokenBonds v_H_outBrokenBonds v_1 dump ${i} all custom 1000 DP_A_DNT_20_${stemp}C id type x y z c_pePerAtom thermo 200 thermo_modify lost warn thermo_modify format 8 %4.0f thermo_modify format 9 %4.0f thermo_modify format 10 %6.0f thermo_modify format 11 %3.0f thermo_modify format 12 %3.0f thermo_modify format 13 %3.0f thermo_modify format 14 %5.0f thermo_modify format 15 %3.0f thermo_modify format 16 %3.0f thermo_modify format 17 %5.0f thermo_modify format 18 %3.0f ##---------------RELAXATION-------------------------------------- if "$i ==1 " then & "run 50000" & else "run 200000" ##---------------BiSection method-------------------------------------- if "${grossBrokenBonds}>=3.333333 && $(abs(v_T2-v_T1)) >1" then & "variable T2 equal ${stemp}" & "variable stemp equal (v_T1+v_T2)/2" & "print 'asigning ${stemp} as the new Temperature & <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<'" & elif "${grossBrokenBonds}<3.333333 && ${T1}>${T2}" & "print 'the system is still stable because the broken bonds =${grossBrokenBonds}'" & "variable T1 equal ${stemp}" & "variable stemp equal ${T1}+400" & elif "${grossBrokenBonds}<3.333333 && ${T1}<${T2} && $(abs(v_T2-v_T1)) >1 && $i>1" & "variable T1 equal ${stemp}" & "variable stemp equal (v_T1+v_T2)/2" & "print 'asigning ${stemp} as the new Temperature & <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<'" & elif "$(abs(v_T2-v_T1)) <1" & "print 'The critical Temperature is = ${stemp} & <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<'" & quit if "$i==1" then & "write_data relaxed_DP_A_DNT_20.data nocoeff nofix" clear next i jump DP_A_DNT_20.in loop ########################################## ##########################################