# Al2O3 crystal, qeq on, minimizes, then calculates elastic constants

variable T_depart equal 300

variable dt equal 0.0002

variable a equal 4.5937
variable c equal 2.9587
variable ca equal ${c}/${a}

variable nx equal 6
variable ny equal 6
variable nz equal 11

variable bx equal ${a}*${nx}
variable by equal ${a}*${ny}
variable bz equal ${c}*${nz}
# =======================================================================

units		metal
atom_style	charge
dimension       3
boundary        p p p


lattice		sc 1.0
region          box_vide prism 0 ${bx} 0 ${by} 0 ${bz} 0.0 0.0 0.0
create_box      2 box_vide

#lattice         sc 1.0
#region          box_TiO2 block 0 ${bx} 0 ${by} 0 ${bz}

# titanium atoms
lattice 	custom ${a} origin 0.0 0.0 0.0 &
                            orient x 1 0 0 orient y 0 1 0 orient z 0 0 1 &
                            a1 1.0  0.0  0.0  &
                            a2 0.0  1.0  0.0  &
                            a3 0.0  0.0 ${ca} &
                            basis 0.0    0.0     0.0 &
 			    basis 0.5    0.5     0.5

create_atoms 	2 region box_vide

# Oxygen atoms
lattice         custom ${a} origin 0.0 0.0 0.0 &
                            orient x 1 0 0 orient y 0 1 0 orient z 0 0 1 &
                            a1 1.0  0.0  0.0  &
                            a2 0.0  1.0  0.0  &
                            a3 0.0  0.0 ${ca} &
                            basis 0.30478    0.30478     0.0 &
                            basis 0.69522    0.69522     0.0 &
			    basis 0.19522    0.80478     0.5 &
			    basis 0.80478    0.19522     0.5 

create_atoms	1 region box_vide
 

mass            1 16.00
group           Oxy type 1
compute         chargeOxy Oxy property/atom q
compute         q_Oxy Oxy reduce ave c_chargeOxy

mass		2 47.867
group    	Ti type 2
compute   	chargeTi Ti property/atom q
compute   	q_Ti Ti reduce ave c_chargeTi

velocity	all create ${T_depart} 277387

pair_style	smtbq
pair_coeff	* * ffield.smtbq.TiO2 O Ti

neighbor	0.5 bin
neigh_modify	every 20 delay 0 check yes

timestep        ${dt}

thermo_style    custom step temp press pe ke etotal c_q_Ti c_q_Oxy lx ly lz vol 
thermo_modify	flush yes
thermo		1


#dump            5 all custom 500 boxAlpha_alumina.lammpstrj id type q x y z  

fix		3 all nve
run		10

unfix		3
#thermo		15
fix             1 all box/relax tri 0.0 vmax 0.001
minimize        1.0e-3 1.0e-5 1000 10000

unfix		1
thermo		1
fix             3 all nve
run             10

