#=========================================================================================
#						(0)							 	                          SETTINGS
#========================================================================================= 
variable			a0		equal 	7.16 
variable			N		equal	2

variable        	Tini	equal   "1750."       # Initial temperature (K)
variable        	Tfin	equal   "5200."       # Melting temperature (K)

variable 			tstep	equal	"0.001" 		# in ps ==> 1 fs
variable 			tdamp 	equal 	"v_tstep*100"  # 
variable 			pdamp 	equal 	"v_tstep*1000" # 
variable 			out1    equal 	"10000"	# Thermo  outputs

#=========================================================================================
#						(I)							                        INITIALIZATION
#=========================================================================================

units           metal	# Units are: A (length), ps (time), eV (nrj), g/cm^3 (density)&
                        # force (eV/atom)
dimension			3
boundary        p p p	# Periodic boundary conditions
atom_style      full	# Affects what quantities are stored for each atom (here only the default values)

#=========================================================================================
#						(II)						                       ATOM DEFINITION
#=========================================================================================
# Creation of a diamond lattice for Si atoms
# It means a fcc structure with four additional atoms inside the cell
# The O atoms are located at average distance between Si atoms
# A diamond lattice can be defined with two basis atoms located at (0,0,0) and a0/4(1,1,1) 
# and three basis vectors a1 a2 a3
# The Si atoms are on the basis atoms (2) and the O in between (6)
# ao, Lattice parameter for beta-cristobalite @ T = 1750 K
# Basis vectors (non normalised) a1 a2 a3
lattice           	custom 	${a0}	& 
					origin	0. 		0. 		0. &
					spacing	1. 		1.		1. &
					a1		0. 		0.5 		0.5 & 
					a2		0.5 		0. 		0.5 &
					a3		0.5 		0.5 		0. &
					basis	0. 		0. 		0. &
					basis 	0.25	0.25	0.25 &
					basis	0.125	0.125	0.125 &
					basis	0.125	0.125	0.625 &
					basis	0.625	0.125	0.125 &
					basis	0.125	0.625	0.125 
 
region				sim_box block 0 ${N} 0 ${N} 0 ${N} units lattice
create_box			2 sim_box 

# Creation of Si atoms on the basis defined in lattice (1&2) and creation of O atoms (3-6)
create_atoms		2 region sim_box &
					basis 1 1 basis 2 1 & 
					basis 3 2 basis 4 2 basis 5 2 basis 6 2 
 
mass				1 28.0855
mass                2 15.9994
group				Si type 1
group				O  type 2
set 				group Si charge 2.4
set					group O  charge -1.2
 
#=========================================================================================
#						(III)						                          FORCE FIELDS
#=========================================================================================

pair_style 		table linear 968
pair_coeff		1 1 pot_SIO2.txt SIO2_1_1 10.17
pair_coeff		1 2 pot_SIO2.txt SIO2_1_2 10.17
pair_coeff		2 2 pot_SIO2.txt SIO2_2_2 10.17

neighbor 2.0 bin
 
neigh_modify delay 0 every 1 check yes

pair_write		1 1 968 r 0.5 10.17 table.txt SiSi
pair_write		1 2 968 r 0.5 10.17 table.txt SiO
pair_write		2 2 968 r 0.5 10.17 table.txt OO

#=========================================================================================
#						(IV)				                          SIMULATION PROCESSES
#=========================================================================================

# Timestep to make the integration of the motion equation (in ps): here 1 fs
timestep        ${tstep}

# Parameters to print out on the screen and log.lammps file
thermo_style    custom  &
				step temp etotal enthalpy vol density press 
thermo          ${out1} # Output every N steps

  
#-----------------------------------------------------------------------------------------
# 						(1) 			                       HEATING FROM 1700 TO 5200 K
#-----------------------------------------------------------------------------------------  
			
# Initial velocities 
velocity        all create ${Tini} 458273253 dist gaussian

# Timestep to make the integration of the motion equation (in ps): here 1 fs
timestep        ${tstep}

# Heating ramp with the NPT ensemble at zero pressure

fix       		1 all npt temp ${Tini} ${Tfin} ${tdamp} iso 0.0 0.0 ${pdamp}
run       		1000000	

write_restart 	restart.heating