# Define unit set and class of atomic model
units metal
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

lattice custom ${initLatConst} a1 1.0 0.0 0.0 a2 -0.5 0.866025403784439 0.0 a3 0.0 0.0 400.0 basis 0.3333333333333333 0.6666666666666667 0.0 basis 0.6666666666666667 0.3333333333333333 0.0 spacing 1.0 0.866025403784439 400.0
region box prism 0 1 0 1 0.0 1.0 -0.5 0 0 units lattice
create_box 1 box 
create_atoms 1 box  
mass * 1.0 

# Specify which KIM Model to use, letting the model compute the virial/pressure
pair_style kim KIMvirial LennardJones612_UniversalShifted__MO_959249795837_002
pair_coeff * * C  

# Set what information to write to dump file
dump id all custom 100 output/lammps.dump id type x y z fx fy fz
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"

# Set what thermodynamic information to print to log
thermo_style custom step atoms lx pe press pxx pyy pxy # 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 "-c_thermo_pe/v_natoms"
variable surften   equal "pxx*lz"
variable a         equal "lx"

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