#  This is the control script for LAMMPS

echo			both

units			metal
boundary		p p p
atom_style		charge


read_data               structure.dat

pair_style              comb3 polar_off
pair_coeff              * * ffield.comb3 Cu O

variable		sysvol		equal vol
variable		sysmass		equal mass(all)/6.0221367e+23
variable		sysdensity	equal v_sysmass/v_sysvol/1.0e-24
variable		etotal		equal etotal
variable		pe		equal pe
variable		ke		equal ke
variable		p		equal press
variable		pxx		equal pxx
variable		pyy		equal pyy
variable		pzz		equal pzz
variable		pyz		equal pyz
variable		pxz		equal pxz
variable		pxy		equal pxy
variable		sxx		equal -pxx
variable		syy		equal -pyy
variable		szz		equal -pzz
variable		syz		equal -pyz
variable		sxz		equal -pxz
variable		sxy		equal -pxy
variable		fmax		equal fmax
variable	        fnorm		equal fnorm
variable		time equal step*dt+0.000001


#-------------------------------------------------------------------------------
# Compute charge vars and build supercell
#-------------------------------------------------------------------------------

replicate 		5 5 5

group 			type1 type 1
compute   		charge1 type1 property/atom q
compute   		q1 type1 reduce ave c_charge1

group    		type2 type 2
compute   		charge2 type2 property/atom q
compute   		q2 type2 reduce ave c_charge2

variable   		qtot equal count(type1)*c_q1+count(type2)*c_q2


#-------------------------------------------------------------------------------
# Set the initial velocities for $T
#-------------------------------------------------------------------------------

velocity		all create 298.15 72489 dist gaussian mom yes rot no

#-------------------------------------------------------------------------------
# NPT 150ps
#-------------------------------------------------------------------------------

reset_timestep		0

thermo_style		custom step v_time press vol v_sysdensity temp pe ke c_q1 c_q2 v_qtot
thermo			100

fix			1 all npt temp 298.15 298.15 0.1 iso 0 0 0.1 drag 0 mtk yes nreset 2000

fix			2 all ave/time 1 999 1000 v_time c_thermo_temp c_thermo_press v_sysvol v_sysdensity v_etotal c_q1 c_q2 v_qtot v_pe v_ke file NPT150ps_averages.txt off 1
fix			3 all ave/time 10  1   10 v_time c_thermo_temp c_thermo_press v_sysvol v_sysdensity v_etotal c_q1 c_q2 v_qtot v_pe v_ke file NPT150ps_instantaneous.txt off 1

fix                     4 all qeq/comb 10 0.001 file fq.out

restart			150000 NPT150ps.restart

dump			trj all custom 1000 NPT150ps.Trajectory.xyz id type q x y z


timestep		0.001

run			150000

undump			trj

unfix			1
unfix			2
unfix			3
unfix			4