# Calculation of properties of alpha quartz using BKS  potential and minimization

units           metal
newton          on
atom_style      charge
dimension       3
boundary        p p p

#lattice custom 5.4054 origin 0.0 0.0 0.0 &
#a1 0.9108 0 0 &
#a2 -0.4554 0.7888 0 &
#
lattice custom 5.4054 origin 0.0 0.0 0.0 &
a1 0.9108 0 0 &
a2 -0.4554 0.7888 0 &
a3 0 0 1 &
basis  0.468       0.0         0.0        &
basis  0.414       0.264       0.121      &
basis  0.0         0.468       0.6667     &
basis  0.532       0.532       0.3333     &
basis  0.736       0.15        0.7877     &
basis  0.85        0.586       0.4543     &
basis  0.264       0.414       0.54562    &
basis  0.15        0.736       0.8789     &
basis  0.586       0.85        0.2123

variable n equal 10
variable mydx equal $n/(1+0.45475/0.9095)

region        simbox block 0 ${mydx} 0 $n 0 $n  units lattice
#region        simbox block 0 $n 0 $n 0 $n  units lattice


create_box 2 simbox
create_atoms 1 region simbox &
basis 1 1 &
basis 2 2 &
basis 3 1 &
basis 4 1 &
basis 5 2 &
basis 6 2 &
basis 7 2 &
basis 8 2 &
basis 9 2


mass 1 28.09
mass 2 16.0

group    siliconatoms type 1
group    oxygenatoms type 2


set    group siliconatoms charge 2.4
set    group oxygenatoms charge -1.2


kspace_style ewald/n 1.0e-5
pair_style buck/coul/long 10.0
pair_coeff 1 1 0.0 0.0001 0.0
pair_coeff 1 2 18003.7572 0.205205 133.5381
pair_coeff 2 2 1338.7730 0.362319 175.0000

# parameters for the thermostat/barostat fix
variable r equal 57085
variable t equal 300 # temperature (in K)
variable p equal 1  # pressure (in bar(s) if units = metal)
variable td equal 0.05 # 100 time steps
variable pd equal 0.5 # 1000 time steps

# initialize
velocity all create $t 28459 rot yes dist gaussian mom yes
reset_timestep 0

dump init all xyz 100000 init_conf_Quartz$n.xyz

# parameters for the thermostat/barostat fix
variable r equal 57085
variable t equal 300 # temperature (in K)
variable p equal 1  # pressure (in bar(s) if units = metal)
variable td equal 0.05 # 100 time steps
variable pd equal 0.5 # 1000 time steps

fix      1 all npt temp $t $t ${td} iso $p $p ${pd}
fix      2 all phonon 10 50000 500000 quartz$nwtf.map.in Quartz$n-nontilted

#
timestep 0.0005

# output
#                    1    2    3  4  5     6   7  8  9  10 11 12
thermo_style  custom step temp pe ke press vol lx ly lz xy xz yz cella cellb cellc cellalpha cellbeta cellgamma
thermo  100000

restart 2000000 restart.one restart.two

# execution
#run 2
run      5500000
write_restart Restart.final