This is the input script (I still have to debug *if command* but regardless of that I am geting error message ERROR: Fix bond/create requires special_bonds lj = 0, 1, 1 (…/fix_bond_create.cpp:187) ) .

units real

boundary p p p

atom_style full

read_restart After_NPT7.restart

pair_style lj/cut/coul/cut 16.0

bond_style harmonic

angle_style harmonic

dihedral_style opls

improper_style harmonic

pair_modify mix geometric tail yes

special_bonds lj/coul 0.0 1.0 1.0 extra 30

pair_coeff 1 1 0.07 3.55

pair_coeff 2 2 0.066 3.5

pair_coeff 3 3 0.066 3.5

pair_coeff 4 4 0.066 3.5

pair_coeff 5 5 0.03 2.42

pair_coeff 6 6 0.03 2.5

pair_coeff 7 7 0 0

pair_coeff 8 8 0.03 2.5

pair_coeff 9 9 0.03 2.5

pair_coeff 10 10 0.0460 0.400

pair_coeff 11 11 0.0005 5.17

pair_coeff 12 12 0.17 3.25

pair_coeff 13 13 0.28 2.93

pair_coeff 14 14 0.25 3.22

pair_coeff 15 15 0.1521 3.1507

pair_coeff 16 16 0.395 3.56

pair_coeff 10 15 0.0836 1.7753

bond_coeff 1 469.00 1.4000

bond_coeff 2 317.00 1.5100

bond_coeff 3 367.00 1.0800

bond_coeff 4 340.00 1.0900

bond_coeff 5 340.00 1.7900

bond_coeff 6 340.00 1.0900

bond_coeff 7 367.00 1.4710

bond_coeff 8 340.00 1.0900

bond_coeff 9 367.00 1.4710

bond_coeff 10 553.00 0.9450

bond_coeff 11 450 0.9572

bond_coeff 12 700.00 1.5300

angle_coeff 1 63.00 120.00

angle_coeff 2 70.00 120.00

angle_coeff 3 35.00 120.00

angle_coeff 4 35.00 109.50

angle_coeff 5 80.00 111.20

angle_coeff 6 62.00 96.00

angle_coeff 7 74.00 107.00

angle_coeff 8 50.00 113.00

angle_coeff 9 50.00 113.00

angle_coeff 10 33.00 107.80

angle_coeff 11 35.00 109.50

angle_coeff 12 33.00 107.80

angle_coeff 13 35.00 109.50

angle_coeff 14 33.00 107.80

angle_coeff 15 35.00 109.50

angle_coeff 16 55 104.52

dihedral_coeff 1 0.000 0.0 7.250 0.0

dihedral_coeff 2 0.000 0.0 7.250 0.0

dihedral_coeff 3 0.000 0.0 7.250 0.0

dihedral_coeff 4 0.000 0.0 0.000 0.0

dihedral_coeff 5 0.000 0.0 0.000 0.0

dihedral_coeff 6 1.438 -0.124 0.264 0.0

dihedral_coeff 7 0.000 0.0 7.250 0.0

dihedral_coeff 8 0.000 0.0 7.250 0.0

dihedral_coeff 9 0.000 0.0 0.350 0.0

dihedral_coeff 10 0.00 0.0 0.350 0.0

dihedral_coeff 11 0.000 0.0 0.302 0.0

dihedral_coeff 12 0.000 0.0 0.302 0.0

dihedral_coeff 13 0.000 0.0 0.302 0.0

dihedral_coeff 14 0.0 5.0 0.0 0.0

dihedral_coeff 15 0.0 5.0 0.0 0.0

neighbor 1.5 bin

neigh_modify every 10 delay 20 check yes

group hydroxide type 7 14

group newbond type 4 14

fix 1 all nve

fix 2 all bond/create 10 4 14 3.5 7

compute 1 newbond property/local btype

compute 2 newbond property/local batom1

compute 3 newbond property/local batom2

run 0

variable a equal c_3-29

variable b equal a+1

variable d equal c_1[1]

compute 6 hydroxide property/atom id

compute 4 hydroxide pe/atom

compute 5 hydroxide ke/atom

run 10000 every 10 “if ‘${d}==0’ then NULL”

" elif ‘(c_4[v_a]+c_4[v_b]c_5[v_a]+c_5[v_b])>110’ " & “delete_bonds reaction atom 4 14 any” &

“elif ‘(c_4[v_a]+c_4[v_b]c_5[v_a]+c_5[v_b])<110’ delete_bonds”& “newbond 4 12”