# Uniaxial deformation tantalum along [100], [110] or [111] # The equilibration step allows the lattice to expand to a temperature of 300 K with a # pressure of 0 kbar at each simulation cell boundary. Then, the simulation cell is deformed in the x-direction at # a strain rate of ???????, while the lateral boundaries are controlled using the NVE equations of motion # to mantain the lenght of the computational box constant. The stress and strain values are output to a separate # file, which can be imported in a graphing application for plotting. The cfg dump files include the x, y, and z # coordinates, the centrosymmetry values, the potential energies, and forces for each atom. This can be directly %# visualized using AtomEye We assume 1 timestep is equal to 0,005 pico-seconds # ---------- Initialize Simulation --------------------- clear units metal dimension 3 boundary p p p atom_style atomic atom_modify map array # ---------- define variables --------------------- variable stemperature equal 300 # temperature in kelvin variable epercentage equal 0.30 # the percentage the body is compressed (CHANGE THIS VALUE) variable myseed equal 12345 # the value seed for the velocity (CHANGE THIS VALUE) variable atomrate equal 250 # the rate in timestep that atoms are dump as CFG (CHANGE THIS VALUE) variable time_step equal 0.005 # time step in pico seconds (CHANGE THIS VALUE) variable tau equal 5 # quasi-isentropic scheme final run time (in time UNITS) variable time_eq equal 500 # time steps for the equlibration part variable tdamp equal "v_time_step*100" # DO NOT CHANGE variable pdamp equal "v_time_step*1000" # DO NOT CHANGE variable R equal "(v_epercentage/0.165)/v_tau" # Compress rate in EMAX (DO NOT CHANGE) variable time_run equal "v_tau/v_time_step" # time step for the deformation part timestep ${time_step} # DO NOT CHANGE # ---------- Create Atoms --------------------- lattice bcc 3.304 region box block 0 10 0 10 0 10 units lattice create_box 1 box # lattice bcc 3.304 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1 lattice bcc 3.304 orient x 1 1 0 orient y 0 0 1 orient z 1 -1 0 # lattice bcc 3.304 orient x 1 1 1 orient y 1 -1 0 orient z 1 1 -2 create_atoms 1 box # ---------- Define Interatomic Potential --------------------- pair_style eam/fs pair_coeff * * Ta1setfl.eam Ta # ---------- Define Settings ---------------------------------- # "compute" Define a computation that will be performed on a group of atoms. Quantities # calculated by a compute are instantaneous values, meaning they are calculated from # information about atoms on the current timestep or iteration, though a compute may internally store some # information about a previous state of the system. Defining a compute does not perform a computation. Instead # computes are invoked by other LAMMPS commands as needed compute myCN all cna/atom bcc compute myKE all ke/atom compute myPE all pe/atom # ---------- Equilibration --------------------------------------- reset_timestep 0 velocity all create ${stemperature} ${myseed} mom yes rot no # The equilibration step NPT allows the lattice to expand to a temperature of 300 K with a # pressure of 0 bar at each simulation cell boundary. The Tdamp and Pdamp parameter is specified in time units # and determines how rapidly the temperature is relaxed. A good choice for many models is a Tdamp of around 100 # timesteps and A good choice for many models is a Pdamp of around 1000 timesteps. Note that this is NOT the same # as 1000 time units for most units settings. fix equilibration all npt temp ${stemperature} ${stemperature} ${tdamp} tchain 1000 iso 0 0 ${pdamp} drag 12 # Output thermodynamics into outputfile # for units metal, pressure is in [bars] = 100 [kPa] = 1/10000 [GPa] # p2, p3, p4 are in GPa variable eq1 equal "step" variable eq2 equal "pxx/10000" variable eq3 equal "pyy/10000" variable eq4 equal "pzz/10000" variable eq5 equal "lx" variable eq6 equal "ly" variable eq7 equal "lz" variable eq8 equal "temp" variable eq9 equal "etotal" fix output1 all print 10 "${eq1} ${eq2} ${eq3} ${eq4} ${eq5} ${eq6} ${eq7} ${eq8} ${eq9}" file run.${stemperature}K.out screen no thermo 10 thermo_style custom step pxx pyy pzz lx ly lz temp etotal # Use cfg for AtomEye (here only print the coordinates of the atoms) dump 1 all cfg ${time_eq} bcc.tantalum.*.cfg id type xs ys zs dump_modify 1 element Ta # RUN AT LEAST 10000 timesteps run ${time_eq} # Store final cell length for strain calculations variable tmp equal "lx" variable L0 equal ${tmp} print "Initial Length, L0: ${L0}" #reset unfix equilibration undump 1 unfix output1 # ---------- Deformation --------------------------------- reset_timestep 0 # In our simulations We seek to control the lateral boundaries Ly and Lz, but using (NEMD) # i.e When the computational box is stretched, the contraction or transverse strain perpendicular to the # load is zero fix 1 all nve fix 2 all deform 1 x qrate ${R} ${tau} units box remap x #IMPORTANT NOTE: When non-equilibrium MD (NEMD) simulations are performed using this fix, the option "remap v" #should normally be used. This is because fix nvt/sllod adjusts the atom positions and velocities to induce a #velocity profile that matches the changing box size/shape. Thus atom coordinates should NOT be remapped by fix #deform, but velocities SHOULD be when atoms cross periodic boundaries, since that is consistent with maintaining #the velocity profile already created by fix nvt/sllod. LAMMPS will warn you if the remap setting is not #consistent with fix nvt/sllod. # Output strain and stress info to file # for units metal, pressure is in [bars] = 100 [kPa] = 1/10000 [GPa] # pxx, pyy, pzz are in GPa variable strain equal "-(lx - v_L0)/v_L0" variable shear equal "0.5*(pxx/10000 - 0.5*(pyy/10000 + pzz/10000))" variable tstep equal "step" variable out1 equal "v_strain" variable out2 equal "pxx/10000" variable out3 equal "pyy/10000" variable out4 equal "pzz/10000" variable out5 equal "v_shear" variable out6 equal "lx" variable out7 equal "ly" variable out8 equal "lz" variable out9 equal "temp" fix def1 all print 1 "${tstep} ${out2} ${out3} ${out4} ${out5} ${out6} ${out7} ${out8} ${out9}" file stress.${stemperature}K.${epercentage}e.out screen no #Use cfg for AtomEye dump 1 all cfg ${atomrate} dump._*.cfg id type xs ys zs c_myPE c_myKE c_myCN fx fy fz dump_modify 1 element Ta # Display thermo variable thermostep equal "v_time_run/10" thermo ${thermostep} thermo_style custom "v_strain" pxx pyy pzz lx ly lz temp etotal pe ke run ${time_run} #${time_run} unfix def1 unfix 1 unfix 2 # ONLY USE THIS STEP ONLY IF YOU KNOW WHERE THE ELASTIC-PLASTIC threshold OCCOURS # FIRST COMPRESS THE CRYSTAL TO THE CRITICAL STRAIN VALUE USING THE PREVIOUS RUNS # (CHANGE VARIABLE) epercentage in #----define variables------# # THE FINAL STEP WILL RUN THE CRYSTAL IN (NVT) AND DUMP .CFG FILES, AND HOPEFULLY THE PLASTIC # BEHAVIOUR WILL APPEAR # -----------FINAL STEP ------------------# # YOU NEED TO UNCOMMENT TO RUN THIS STEP # here we run in the NVT enssamble and we look for elastic-pleastic transition. # that's easy to implement # fix 1 all nvt temp ${stemperature} ${stemperature} 0.1 drag 1 # variable p1 equal "step" # variable p2 equal "-pxx/10000" # variable p3 equal "-pyy/10000" # variable p4 equal "-pzz/10000" # variable p5 equal "lx" # variable p6 equal "ly" # variable p7 equal "lz" # variable p8 equal "temp" # variable p9 equal "etotal" # fix output1 all print 10 "${p1} ${p2} ${p3} ${p4} ${p5} ${p6} ${p7} ${p8} ${p9}" file nvt.${stemperature}K.${epercentage}e.out screen no # thermo 1000 # thermo_style custom step pxx pyy pzz lx ly lz temp etotal # dump 1 all cfg 250 dump._*.cfg id type xs ys zs c_peratom c_cnaatom fx fy fz # dump_modify 1 element Ta # run 20000 # run by 20 ps # SIMULATION DONE clear print "creo ya esta =)"