# ------------------------ INITIALIZATION ---------------------------- units real # real = grams/mole, angstroms, femtosecs, Kcal/mole, Kelvin, atmospheres dimension 3 # 3D boundary p p p # periodic in x,y,z atom_style charge # atomic = only default values, e.g. no charge, … # ---------- Define Variables ---------------- shell mkdir output variable latparam equal 2.846 #2.86 variable repeat equal 1 # ----------------------- ATOM DEFINITION ---------------------------- lattice bcc ${latparam} region whole block 0 10 0 10 0 10 create_box 1 whole lattice bcc ${latparam} orient x 1 0 0 orient y 0 1 0 orient z 0 0 1 create_atoms 1 region whole replicate ${repeat} ${repeat} ${repeat} # Replicate domain x,y,z mass 1 55.8450 # ------------------------ FORCE FIELDS ------------------------------ pair_style reax/c NULL # Set the formula(s) LAMMPS uses to compute pairwise interactions pair_coeff * * FeCrOSCH.reaxff Fe # Specify the pairwise force field coefficients for one or more pairs of atom types; * * means for all atom types neighbor 2.0 bin neigh_modify delay 0 every 10 check yes # ------------------------- Define Computes --------------------------------- compute csym all centro/atom bcc # centrosymmetry parameter CS = SUM(i=1 to N/2)[|R_i + R_{i+N/2}|^2] compute eng all pe/atom # c_eng defined over all atoms to compute potential energy per atom ###################################### # EQUILIBRATION # reset timer reset_timestep 0 # 2 fs time step timestep 2 # initial velocities velocity all create 300 12345 mom yes rot no # temperature random_seed # If mom = yes, linear momentum of the newly created ensemble of velocities is zeroed; # If rot = yes, angular momentum is zeroed. # thermostat + barostat fix 1 all npt temp 300.0 300.0 1000.0 iso 0.0 0.0 1000.0 drag 1.0 fix 2 all qeq/reax 1 0.0 10.0 1e-6 param.qeq # instrumentation and output variable s1 equal "time" variable s2 equal "lx" variable s3 equal "ly" variable s4 equal "lz" variable s5 equal "vol" variable s6 equal "press" variable s7 equal "pe" variable s8 equal "ke" variable s9 equal "etotal" variable s10 equal "temp" fix writer all print 250 "${s1} ${s2} ${s3} ${s4} ${s5} ${s6} ${s7} ${s8} ${s9} ${s10}" file output/Fe_eq.txt screen no # thermo thermo 50 thermo_style custom step time cpu cpuremain lx ly lz press pe temp # dumping trajectory dump 1 all atom 250 output/Fe_bcc_tensile_eq.dump # 24 ps MD simulation (assuming 2 fs time step) run 12000 # clearing fixes and dumps unfix 1 unfix 2 undump 1 # saving equilibrium length for strain calculation variable tmp equal "lx" variable L0 equal ${tmp} print "Initial Length, L0: ${L0}" ###################################### # DEFORMATION # reset timer reset_timestep 0 # 2 fs time step timestep 2 # thermostat + barostat # temp Tstart Tstop Tdamp. The Tdamp parameter is specified in time units and determines how rapidly the temperature is relaxed. # For example, a value of 10.0 means to relax the temperature in a timespan of (roughly) 10 time units fix 1 all npt temp 300.0 300.0 1000.0 y 0 0 1 z 0 0 1 drag 1.0 # fix 2 all qeq/reax 1 0.0 10.0 1e-6 param.qeq # nonequilibrium straining in x-direction at strain rate = 1x10^10 / s = 1x10^-5 / fs in units real variable srate equal 1.0e10 variable srate1 equal "v_srate / 1.0e15" fix 3 all deform 1 x erate ${srate1} units box remap x # deform N x = perform box deformation every N timesteps in x direction # erate = engineering strain rate (1/time units) # instrumentation and output # for units real, pressure is in [atmospheres] = 101.325 [kPa] = 0.000101325 [GPa] => p2, p3, p4 are in GPa variable strain equal "(lx - v_L0)/v_L0" variable p1 equal "v_strain" variable p2 equal "-pxx*0.000101325" variable p3 equal "-pyy*0.000101325" variable p4 equal "-pzz*0.000101325" fix writer all print 125 "${p1} ${p2} ${p3} ${p4}" file output/Fe_deform.txt screen no # thermo thermo 1 thermo_style custom step cpuremain v_strain v_p2 v_p3 v_p4 press pe temp # dumping standard atom trajectories dump 1 all atom 125 output/Fe_bcc_tensile_deform.dump # dumping custom cfg files containing coords + ancillary variables dump 2 all cfg 125 output/Fe_bcc_tensile_deform_*.cfg mass type xs ys zs c_csym c_eng fx fy fz dump_modify 2 element Fe # 20 ps MD simulation (assuming 2 fs time step) run 10000 # clearing fixes and dumps unfix 1 unfix 2 unfix 3 unfix writer undump 1 undump 2 ###################################### # SIMULATION DONE print "All done!"