#================================ Simulation Setting ================================ # input data file variable lmp_data string data.alni # control file variable lmp_control string input_lmp_control # temperature(K) variable T equal "1728K" #variable T1 equal "2500.00" #variable T2 equal "2500.00" # timestep(fs) variable dt1 equal "0.1000" #variable dt2 equal "0.2000" # number of steps variable nstep1 equal "3000000" #variable nstep2 equal "500000" # thermo output every variable freq_thermo1 equal "1000" #variable freq_thermo2 equal "500" # frame output every variable freq_frame1 equal "1000" #variable freq_frame2 equal "500" # xyz output every for safety (rerun) variable freq_xyz1 equal "1000" #variable freq_xyz2 equal "500" #==================================================================================== units real boundary p p p #s means shrink-wrapped atom_style charge read_data ${lmp_data} pair_style reax/c NULL checkqeq yes pair_coeff * * WaterCaSiAlS.ff.txt Al fix species all reax/c/species 1 1 ${freq_frame1} out.species element Al fix fqeq all qeq/reax 1 0.0 10.0 1e-6 reax/c compute reax all pair reax/c variable eb equal c_reax[1] # bond energy variable ea equal c_reax[2] # atom energy variable elp equal c_reax[3] # lone-pair energy variable emol equal c_reax[4] # molecule energy (always 0.0) variable ev equal c_reax[5] # valence angle energy variable epen equal c_reax[6] # double-bond valence angle penalty variable ecoa equal c_reax[7] # valence angle conjugation energy variable ehb equal c_reax[8] # hydrogen bond energy variable et equal c_reax[9] # torsion energy variable eco equal c_reax[10] # conjugation energy variable ew equal c_reax[11] # van der Waals energy variable ep equal c_reax[12] # Coulomb energy variable efi equal c_reax[13] # electric field energy (always 0.0) variable eqeq equal c_reax[14] # charge equilibration energy neighbor 2 bin neigh_modify every 10 delay 0 check no #neigh_modify delay 0 every 1 check yes page 10000 one 100 #balance 1.0 shift xy 100 1.0 out tmp.balance_xy #balance 1.0 shift xz 100 1.0 out tmp.balance_xz #balance 1.0 shift yz 100 1.0 out tmp.balance_yz #============================== Energy minimization =============================== #thermo 1 #thermo_style custom step cpuremain temp epair etotal enthalpy press & # v_eb v_ea v_elp v_emol v_ev v_epen v_ecoa & # v_ehb v_et v_eco v_ew v_ep v_efi v_eqeq #fix boxrelax all box/relax aniso 0.0 vmax 0.001 #min_style cg #minimize 1.0e-6 1.0e-8 10000 100000 #unfix boxrelax #================================= fix_NVT simulation ================================== #NVT simulation with Nose-Hoover thermostat fix fnptalni all npt temp 1500 1500 100.0 iso 0.0 0.0 1000.0 #NVT simulation with Berendsen thermostat #group nm id 1:1757 #group updown id 1758:4093 #fix fixcnt updown move linear 0 0 0 #fix fnvtB11 nm nve #fix fnvtB12 nm temp/berendsen 2500 2500 100.0 #==================================================================================== thermo ${freq_thermo1} thermo_style custom step cpuremain temp epair etotal enthalpy press & v_eb v_ea v_elp v_emol v_ev v_epen v_ecoa & v_ehb v_et v_eco v_ew v_ep v_efi v_eqeq dump dtrj all custom ${freq_frame1} ${lmp_data}.lammpstrj_npt id type x y z q dump dxyz all xyz ${freq_xyz1} ${lmp_data}_npt_*.xyz dump dsnap all custom ${freq_xyz1} ${lmp_data}_npt_*.lammpstrj id type x y z q dump_modify dxyz pad 10 element Al reset_timestep 0 timestep ${dt1} run ${nstep1} unfix fnvtNH undump dtrj undump dxyz undump dsnap write_restart restart.reaxff unfix fqeq print "all done"