# The units we assume we'll use. The kim-lammps-preprocessor may swap this line out if
# running against a Simulator Model whose units are not 'metal'
units metal

# Variables that can be adjusted by kim-lammps-preprocessor to switch unit sets for
# Simulator Models
variable _u_distance equal 1.0
variable _u_energy equal 1.0
variable _u_force equal 1.0
variable _u_pressure equal 1.0

# Isolated atom energy for this species in eV (computed in a separate LAMMPS calculation)
variable isolated_atom_energy_sum equal -0.195353

# This line may be swapped out by kim-lammps-preprocessor if running against a Simulator
# Model whose atom_style is not 'atomic'
atom_style atomic

# Periodic boundary conditions, z direction big enough so images don't interact
boundary p p p

# Create a hexagonal lattice with the lattice spacing and basis atoms supplied by the
# user using a single primitive hexagonal unit cell. Insert a graphene-like or a TMD
# basis.
variable initLatConst equal 2.5*${_u_distance}

lattice custom ${initLatConst} a1       1.0                0.0    0.0 &
                               a2      -0.5  0.866025403784439    0.0 &
                               a3       0.0                0.0  150.0 &
                               spacing  1.0  0.866025403784439  150.0 &
                               basis 0.3333333333333333 0.6666666666666667 0.0 basis 0.6666666666666667 0.3333333333333333 0.0

region box prism 0 1 0 1 0 1 -0.5 0 0 units lattice
create_box 1 box
create_atoms 1 box  

# Mass unimportant since we're not doing time integration
mass * 1.0

# Specify which KIM Model to use
pair_style kim Sim_LAMMPS_ReaxFF_SinghSrinivasanNeekAmal_2013_CFH__SM_306840588959_000
pair_coeff * * C  

# Variables used to rescale the positions and forces so that the quantities in the
# dumpfile are in the original metal units (angstrom and eV/angstrom) even if we're
# running with a Simulator Model that uses different units
variable  x_metal   atom  x/${_u_distance}
variable  y_metal   atom  y/${_u_distance}
variable  z_metal   atom  z/${_u_distance}
variable  fx_metal  atom    fx/${_u_force}
variable  fy_metal  atom    fy/${_u_force}
variable  fz_metal  atom    fz/${_u_force}

# Set what information to write to dump file
dump id all custom 100 output/lammps.dump id type v_x_metal v_y_metal v_z_metal &
    v_fx_metal v_fy_metal v_fz_metal
dump_modify id format line "%d %d %16.7f %16.7f %16.7f %16.7f %16.7f %16.7f"

# Set boundary conditions to be stress-free
fix 1 all box/relax x 0.0 y 0.0 z 0.0 couple xy

# fix atoms
fix 2 all setforce 0.0 0.0 0.0

# For TMDs, need difference in z-position of layers
#    group  transmetal        type  1
#    group   chalcogen        type  2
#  compute   upperZpos   chalcogen  reduce  max  z
#  compute  middleZpos  transmetal  reduce  max  z

# Consistent with the property definition, we are defining the z position in terms
# of the 2D lattice constant
# variable zpos equal "(c_upperZpos-c_middleZpos)/lx"

# Variables used to rescale the energy and stress so that the quantities in the thermo
# output are in the original metal units (eV and bars) even if we're running with a
# Simulator Model that uses different units
variable     pe_metal  equal      "c_thermo_pe/v__u_energy"
variable     lx_metal  equal              lx/${_u_distance}
variable     lz_metal  equal              lz/${_u_distance}
variable  press_metal  equal "c_thermo_press/v__u_pressure"
variable    pxx_metal  equal             pxx/${_u_pressure}
variable    pyy_metal  equal             pyy/${_u_pressure}
variable    pxy_metal  equal             pxy/${_u_pressure}

# Set what thermodynamic information to print to log
thermo_style custom step atoms v_pe_metal v_lx_metal v_lz_metal &
    press v_press_metal v_pxx_metal v_pyy_metal v_pxy_metal &
    # v_zpos

thermo 100 # Print every 100 timesteps

# Perform static minimization using the Polack-Ribiere conjugate gradient method.
# Restart minimization after end with smaller allowed volume step size
minimize 1e-12 1e-12 10000 50000

# Define auxiliary variables to contain cohesive energy and equilibrium lattice constant
variable    natoms  equal  "count(all)"
variable ecohesive  equal  "-(v_pe_metal - v_isolated_atom_energy_sum)/v_natoms"
variable   surften  equal  "v_pxx_metal*v_lz_metal"
variable         a  equal  "v_lx_metal"

# Output cohesive energy and equilibrium lattice constant
# print "Free z parameter in layer group = ${zpos} angstrom"
print "Final surface tension = ${surften} bar*angstrom"
print "Cohesive energy = ${ecohesive} eV/atom"
print "Equilibrium lattice constant = ${a} angstrom"