# This is the control script for LAMMPS. Run LAMMPS on GPU
echo both
log 2.1_Initialize.out
#-------------------------------------------------------------------------------
# Stage 2.1: Initialize LAMMPS run for 3-d periodic
#-------------------------------------------------------------------------------
units real
boundary p p p
atom_style full
atom_modify map array
box tilt large
read_data structure.dat
pair_style reax/c NULL checkqeq no safezone 10.0 mincap 1000
pair_coeff * * parameters.dat H O Si
variable R equal 0.00198722
variable sysvol equal vol
variable sysmass equal mass(all)/6.0221367e+23
variable sysdensity equal v_sysmass/v_sysvol/1.0e-24
variable coulomb equal ecoul+elong
variable etotal equal etotal
variable pe equal pe
variable ke equal ke
variable evdwl equal evdwl
variable epair equal epair
variable ebond equal ebond
variable eangle equal eangle
variable edihed equal edihed
variable eimp equal eimp
variable lx equal lx
variable ly equal ly
variable lz equal lz
variable Nthermo equal 0
variable cella equal lx
variable cellb equal sqrt(ly*ly+xy*xy)
variable cellc equal sqrt(lz*lz+xz*xz+yz*yz)
variable cellalpha equal acos((xy*xz+ly*yz)/(v_cellb*v_cellc))
variable cellbeta equal acos(xz/v_cellc)
variable cellgamma equal acos(xy/v_cellb)
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
variable surfacetension equal 0.5*v_lz*(0.5*(v_sxx+v_syy)-v_szz)
thermo_style custom step v_time press vol v_sysdensity temp ebond eangle edihed eimp evdwl ecoul etail elong pe ke
thermo_modify flush yes
#
# Set up the fixed and movable groups
#
group movable union all
group fixed subtract all movable
compute sum_f1 movable reduce sum fx fy fz
variable sum_fx equal c_sum_f1[1]
variable sum_fy equal c_sum_f1[2]
variable sum_fz equal c_sum_f1[3]
#
# Subsets
#
group subset_Surface_O id 80 92 104 188 200 212 296 308 320 404
group subset_Surface_O id 416 428 512 524 536 620 632 644 728 740
group subset_Surface_O id 752 836 848 860 944 956 968 1052 1064 1076
group subset_Surface_O id 1160 1172 1184 1268 1280 1292 1376 1388 1400 1484
group subset_Surface_O id 1496 1508 1592 1604 1616 1700 1712 1724 1808 1820
group subset_Surface_O id 1832 1916 1928 1940 2024 2036 2048 2132 2144 2156
group subset_Surface_O id 2240 2252 2264 2348 2360 2372 2456 2468 2480 2564
group subset_Surface_O id 2576 2588 2672 2684 2696 2780 2792 2804 2888 2900
group subset_Surface_O id 2912 2996 3008 3020 3104 3116 3128 3212 3224 3236
group subset_Surface_O id 3320 3332 3344 3428 3440 3452 3536 3548 3560 3644
group subset_Surface_O id 3656 3668 3752 3764 3776 3860 3872 3884 3968 3980
group subset_Surface_O id 3992 4076 4088 4100 4184 4196 4208 4292 4304 4316
group subset_Surface_O id 4400 4412 4424 4508 4520 4532 4616 4628 4640 4724
group subset_Surface_O id 4736 4748 4832 4844 4856 4940 4952 4964 5048 5060
group subset_Surface_O id 5072 5156 5168 5180 5264 5276 5288 5372 5384 5396
log 2.2_NVT.out
#-------------------------------------------------------------------------------
# Stage 2.2: NVT integration for 1000 ps with a timestep of 0.2 fs
# Temperature 300 K
#-------------------------------------------------------------------------------
reset_timestep 0
thermo_style custom step v_time press vol v_sysdensity temp ebond eangle edihed eimp evdwl ecoul etail elong pe ke
thermo ${Nthermo}
fix 1 movable nvt temp 300 300 100 drag 0.0
fix 2 movable ave/time 1 500000 500000 v_time c_thermo_temp c_thermo_press v_sysvol v_sysdensity v_etotal v_pe v_ke v_evdwl v_coulomb v_sxx v_syy v_szz v_syz v_sxz v_sxy file 2.2_averages.txt off 1
fix 3 movable ave/time 500 1 500 v_time c_thermo_temp c_thermo_press v_sysvol v_sysdensity v_etotal v_pe v_ke v_evdwl v_coulomb v_sxx v_syy v_szz v_syz v_sxz v_sxy file 2.2_instantaneous.txt
fix 4 all qeq/reax 1 0.0 10.0 1.0e-6 reax/c
restart 5000000 2.2.restart
dump trj all custom 50000 2.2.Trajectory.xyz id mol type q xs ys zs
fix trjE all ave/time 50000 1 50000 v_etotal v_pe v_ke file 2.2.energies.txt
dump sci all custom 500000 2.2.xyz id mol type q xs ys zs
timestep 0.2
run 5000000
undump trj
unfix trjE
undump sci
restart 0
dump sci all custom 5000000 2.2.xyz id mol type q xs ys zs
run 0
undump sci
unfix 1
unfix 2
unfix 3
unfix 4