#-------------define simulation box-------------#
units   metal
atom_style molecular

boundary p p p

#-------------define force dield----------------#

pair_style lj/cut 14.550
bond_style hybrid harmonic morse
angle_style hybrid harmonic cosine/squared
dihedral_style hybrid harmonic multi/harmonic

read_restart restart.polymer.2800000
pair_coeff 1 1 0.00491188 4.28 15
pair_coeff 1 2 0.00491188 4.28 15
pair_coeff 1 3 0.013216 3.87 15
pair_coeff 1 4 0.013216 3.87 15
pair_coeff 1 5 0.013216 3.87 15
pair_coeff 1 6 0.013216 3.87 15
pair_coeff 1 7 0.013216 3.87 15
pair_coeff 1 8 0.013216 3.87 15
pair_coeff 1 9 0.013216 3.87 15
pair_coeff 1 10 0.013216 3.87 15
pair_coeff 2 2 0.00491188 4.28 15
pair_coeff 2 3 0.013216 3.87 15
pair_coeff 2 4 0.013216 3.87 15
pair_coeff 2 5 0.013216 3.87 15
pair_coeff 2 6 0.013216 3.87 15
pair_coeff 2 7 0.013216 3.87 15
pair_coeff 2 8 0.013216 3.87 15
pair_coeff 2 9 0.013216 3.87 15
pair_coeff 2 10 0.013216 3.87 15
pair_coeff 3 3 0 3.46 15
pair_coeff 3 4 0.035559 3.46 15
pair_coeff 3 5 0.035559 3.46 15
pair_coeff 3 6 0.035559 3.46 15
pair_coeff 3 7 0.035559 3.46 15
pair_coeff 3 8 0.035559 3.46 15
pair_coeff 3 9 0.035559 3.46 15
pair_coeff 3 10 0.035559 3.46 15
pair_coeff 4 4 0 3.46 15
pair_coeff 4 5 0.035559 3.46 15
pair_coeff 4 6 0.035559 3.46 15
pair_coeff 4 7 0.035559 3.46 15
pair_coeff 4 8 0.035559 3.46 15
pair_coeff 4 9 0.035559 3.46 15
pair_coeff 4 10 0.035559 3.46 15
pair_coeff 5 5 0 3.46 15
pair_coeff 5 6 0.035559 3.46 15
pair_coeff 5 7 0.035559 3.46 15
pair_coeff 5 8 0.035559 3.46 15
pair_coeff 5 9 0.035559 3.46 15
pair_coeff 5 10 0.035559 3.46 15
pair_coeff 6 6 0 3.46 15
pair_coeff 6 7 0.035559 3.46 15
pair_coeff 6 8 0.035559 3.46 15
pair_coeff 6 9 0.035559 3.46 15
pair_coeff 6 10 0.035559 3.46 15
pair_coeff 7 7 0 3.46 15
pair_coeff 7 8 0.035559 3.46 15
pair_coeff 7 9 0.035559 3.46 15
pair_coeff 7 10 0.035559 3.46 15
pair_coeff 8 8 0 3.46 15
pair_coeff 8 9 0.035559 3.46 15
pair_coeff 8 10 0.035559 3.46 15
pair_coeff 9 9 0 3.46 15
pair_coeff 9 10 0.035559 3.46 15
pair_coeff 10 10 0 3.46 15

mass 1 16
mass 2 16
mass 3* 56

bond_coeff     1 harmonic 20.729 1.53
bond_coeff     2 morse 8.515845 1.55 2.8
bond_coeff     3 harmonic 20.729 3.7
angle_coeff    1 cosine/squared 5.389  112.813
angle_coeff    2 harmonic 17.75327 120
dihedral_coeff 1 multi/harmonic 0.0915373 0.1874586 0.050578 -0.32958 0
dihedral_coeff 2 harmonic 0.17996  -1 2
special_bonds  lj 0 0 0 angle yes dihedral yes


neighbor     2 bin
neigh_modify   delay 10

#---------------define output parameters---------------#

compute         t_pe  all pe 
compute         p_pe  all pe pair
compute         b_pe  all pe bond
compute         a_pe  all pe angle
compute         d_pe  all pe dihedral
compute         t_pea all pe/atom 
compute         p_pea all pe/atom pair
compute         b_pea all pe/atom bond
compute         a_pea all pe/atom angle
compute         d_pea all pe/atom dihedral
compute         astress all stress/atom NULL
compute         m_bond all bond/local dist 
compute         n_bond all property/local batom1 batom2 btype
compute         m_angle all angle/local theta
compute         n_angle all property/local aatom1 aatom2 aatom3 atype
compute         m_dihedral all dihedral/local phi
compute         n_dihedral all property/local datom1 datom2 datom3 datom4 dtype

variable        n equal count(all)
variable        p4 equal -pxx*0.1
variable        p5 equal -pyy*0.1
variable        p6 equal -pzz*0.1
variable        p7 equal c_t_pe
variable        p8 equal c_p_pe
variable        p9 equal c_b_pe
variable        p10 equal c_a_pe
variable        p11 equal c_d_pe
variable        p12 equal -press*0.1
variable        sigma_xx atom c_astress[1]*$n/(vol*10)
variable        sigma_yy atom c_astress[2]*$n/(vol*10)
variable        sigma_zz atom c_astress[3]*$n/(vol*10)
variable        sigma_xy atom c_astress[4]*$n/(vol*10)
variable        sigma_xz atom c_astress[5]*$n/(vol*10)
variable        sigma_yz atom c_astress[6]*$n/(vol*10)
variable        sigma_von atom sqrt(0.5)*sqrt((c_astress[1]-c_astress[2])^2+(c_astress[2]-c_astress[3])^2+(c_astress[3]-c_astress[1])^2+6*(c_astress[4]^2+c_astress[5]^2+c_astress[6]^2))*$n/(vol*10)
variable        st  equal "step"

#-----------------energy minimization------------#

dump         2 all custom 1000 polymer_min.* id mol type mass x y z c_t_pea c_p_pea c_b_pea c_a_pea c_d_pea
dump_modify  2 sort id
dump_modify  2 first yes
fix          min all print 1000 "${st} ${p7} ${p8} ${p9} ${p10} ${p11}" file min_polymer.txt screen no
minimize     1e-12 1e-12 20000 20000
undump       2
unfix        min
#---------------define strain---------------

variable        tmp equal "lx"

variable        L0 equal ${tmp}

variable        strain equal "lx-v_L0"


variable        tmp_y equal "ly"

variable        L0_y equal ${tmp_y}

variable        strain_y equal "ly-v_L0_y"


variable        tmp_z equal "lz"

variable        L0_z equal ${tmp_z}

variable        strain_z equal "lz-v_L0_z"


variable        p1 equal "v_strain/v_L0"

variable        p2 equal "v_strain_y/v_L0_y"

variable        p3 equal "v_strain_z/v_L0_z"

#-------------------heat the sample-----------------------#

reset_timestep  0
dump            1 all local 5000 bond.* index c_n_bond[3] c_m_bond c_n_bond[1] c_n_bond[2]
dump_modify     1 first yes
dump            2 all local 5000 angle.* index c_n_angle[4] c_m_angle c_n_angle[1] c_n_angle[2] c_n_angle[3] 
dump_modify     2 first yes
dump            3 all local 5000 dihedral.* index c_n_dihedral[5] c_m_dihedral c_n_dihedral[1] c_n_dihedral[2] c_n_dihedral[3] c_n_dihedral[4]
dump_modify     3 first yes
dump            4 all custom 5000 polymer_equ.* id mol type mass x y z c_t_pea c_p_pea c_b_pea c_a_pea c_d_pea v_sigma_xx v_sigma_yy v_sigma_zz v_sigma_xy v_sigma_xz v_sigma_yz
dump_modify     4 sort id
dump_modify     4 first yes
fix             equ all print 1000 "${st} ${p7} ${p8} ${p9} ${p10} ${p11} ${p12}" file equ_polymer.txt screen no

timestep        0.001
thermo          100

fix             break_cl all bond/break 1000 3 3.8
fix             1 all npt temp 200.0 200.0 0.1 x 1.0 1.0 1 y 1.0 1.0 1 z 1.0 1.0 1
run             400000
unfix           1
undump          1

undump          2

undump          3

undump          4
#---------------define strain---------------

variable        tmp equal "lx"

variable        L0 equal ${tmp}

variable        strain equal "lx-v_L0"


variable        tmp_y equal "ly"

variable        L0_y equal ${tmp_y}

variable        strain_y equal "ly-v_L0_y"


variable        tmp_z equal "lz"

variable        L0_z equal ${tmp_z}

variable        strain_z equal "lz-v_L0_z"


variable        p1 equal "v_strain/v_L0"

variable        p2 equal "v_strain_y/v_L0_y"

variable        p3 equal "v_strain_z/v_L0_z"

#---------------tensile test---------------
variable         m_step equal 10000
reset_timestep   0
fix              1 all npt temp 200.0 200.0 0.1 y 1 1 1 z 1 1 1

fix              def1  all ave/time 1 500 ${m_step} v_st v_p1 v_p4 v_p5 v_p6 v_p7 v_p8 v_p9 v_p10 v_p11 file tension_polymer.txt off 2 #every 0.001
thermo           100

dump             2 all  custom  ${m_step} polymer_ten.* id mol type mass x y z vx vy vz v_sigma_xx v_sigma_yy v_sigma_zz v_sigma_xy v_sigma_xz v_sigma_yz v_sigma_von c_t_pea c_p_pea c_b_pea c_a_pea c_d_pea
dump_modify      2 sort id   
dump_modify      2 first yes
label            tensile
variable         i loop 1 100
label            dump_or_not
variable         j loop 1 6
variable         x_f equal ${L0}*2/500

fix              stretch all deform ${m_step} x delta 0 ${x_f} remap x units box
run              ${m_step}
unfix            stretch
if               $j==5 then "jump in.polymer break"
next             j

jump             in.polymer dump_or_not
label            break
variable         j delete
write_data       strain_${st}.data
next             i
jump             in.polymer tensile