# ------------------------ Initialization ----------------------------
units       lj

dimension       2

boundary        p   p   p

atom_style      atomic

atom_modify map yes

# ------------------------Read dump------------------------------
read_data       @data.cooling

# ------------------------ Force field------------------------------
variable    sigmall     equal       1.0
variable    sigmass     equal       "5.0/6.0*v_sigmall"
variable    sigmals     equal       "(v_sigmall+v_sigmass)/2.0"

variable    rc        equal       1.5
variable    epsilon       equal       1.0

pair_style      table    linear     750
pair_coeff      1   1       ./Potential/LJll.table  INT ${rc}
pair_coeff      1   2       ./Potential/LJls.table  INT ${rc}
pair_coeff      2   2       ./Potential/LJss.table  INT ${rc}

mass        1       2.0                #Large
mass        2       1.0                 #Small

#-------------------------------unit conversion -------------------------------
variable        sigma       equal       "2.7*10^(-10)"      #m
variable        epsilon     equal       "0.151*1.602176565*10^(-19)"    #J
variable        mass        equal       "46*1.66053886*10^(-27)"    #kg

variable        hightempstar        equal       "2000.0*1.380649*10^(-23)/v_epsilon"    #2000K
variable        lowtempstar         equal       "50.0*1.380649*10^(-23)/v_epsilon"  #50K
variable        pressurestar        equal       "1.0*10.0^9*(v_sigma)^3/v_epsilon"    #5 GPa
variable        timestar        equal       "0.001*10^(-12)*(v_epsilon/v_mass/v_sigma^2)^(1/2)" #0.001ps

# ------------------------ Default setting -------------------------------
timestep        1e-4

thermo      10000
thermo_style        custom      step temp vol pe ke press enthalpy epair

fix     2d       all     enforce2d

# ------------------------ Data Output ------------------------------
variable        time        equal       "time"
variable        step        equal       "step"
variable        temp        equal       "temp"
variable        vol         equal       "vol"
variable        pe          equal       "pe"
variable        ke          equal       "ke"
variable        press       equal       "press"
variable        enthalpy    equal       "enthalpy"
variable        epair       equal       "epair"

#   Deformed box length
variable        LX  equal   "lx"
variable        LY  equal   "ly"
variable        LZ  equal   "lz"

#   Initial box length
variable        LX0     equal       ${LX}
variable        LY0     equal       ${LY}
variable        LZ0     equal       ${LZ}

#   Engineering Strain
variable        strainX     equal       "(lx-v_LX0)/v_LX0"
variable        strainY     equal       "(ly-v_LY0)/v_LY0"
variable        strainZ     equal       "(lz-v_LZ0)/v_LZ0"
variable        strainXY        equal       "xy/v_LY0"
variable        strainXZ        equal       "xz/v_LZ0"
variable        strainYZ        equal       "yz/v_LZ0"

#   Stress unit=GPa
variable        pXX     equal       "-pxx*v_epsilon/(v_sigma)^3/10^9"
variable        pYY     equal       "-pyy*v_epsilon/(v_sigma)^3/10^9"
variable        pZZ     equal       "-pzz*v_epsilon/(v_sigma)^3/10^9"
variable        pXY     equal       "-pxy*v_epsilon/(v_sigma)^3/10^9"
variable        pXZ     equal       "-pxz*v_epsilon/(v_sigma)^3/10^9"
variable        pYZ     equal       "-pyz*v_epsilon/(v_sigma)^3/10^9"

#	Atomic stress
compute       vor     all     voronoi/atom
compute       stress      all     stress/atom     thermo_temp

variable      Sxx     atom        c_stress[1]/c_vor[1]*v_epsilon/(v_sigma)^3/10^9
variable      Syy     atom        c_stress[2]/c_vor[1]*v_epsilon/(v_sigma)^3/10^9
# variable      Szz     atom        c_stress[3]/c_vor[1]*v_epsilon/(v_sigma)^3/10^9
variable      Sxy     atom        c_stress[4]/c_vor[1]*v_epsilon/(v_sigma)^3/10^9
# variable      Sm      atom        (v_Sxx+v_Syy+v_Szz)/3.0
variable      Sm      atom        (v_Sxx+v_Syy)/2.0

variable    n_tol  equal   200
variable    n       loop    ${n_tol}

fix     data        all print       1        "${n}	${time}	${strainX}	${strainY}	${strainZ}	${strainXY}	${strainXZ}	${strainYZ}	${pXX}	${pYY}	${pZZ}	${pXY}	${pXZ}	${pYZ}"     file        @shear.txt     screen no

# ------------------------ Ensemble -----------------------------

# fix     ensemble        all     nvt     temp    ${lowtempstar}  ${lowtempstar}  $(100.0*dt)

# ------------------------ Loading Definition -----------------------------

change_box      all     triclinic

variable        strain_rate     equal       2
### ------------------------ Loop -----------------------------

label loop

variable    a   equal   $n%1

fix     loading     all     deform      1   xy      erate   ${strain_rate}    units   box     remap   x       flip    yes

run     1

unfix       loading

# ------------------------ Energy minimization -----------------------------

if "$a == 0" then "dump 1 all custom 10000000 load/dump.$n id type x y z" &
"dynamical_matrix all regular 0.01 file matrix/dynmat.$n.dat binary no"

minimize 0.0 1.0e-8 100000 10000000

# if "$a == 0" then "write_data  data.minimize.$n"

if "$a == 0" then "undump 1"

next n

jump @AQS.in loop

# ------------------------ Done ------------------------------
print       "All done!"