I have modeled the bending of a silicon cantilever nanobeam in the lammps, using the commands addforce (fix 3 Gr_Bending_strip addforce 0 ${loadforce} 0)
The code works.
But my goal is to plot the curve of applied force against the deflection of its free end and also to obtain Young’s modulus and I do not know how to do that. Could anybody help me and instruct me on how to obtain these data? How can I calculate the applied force?
This is the code I wrote:
###INPUTS###
variable dt equal 0.001
variable loadforce equal -0.006
variable T equal 300
variable iterequi equal 00000
variable iterrun equal 700000
Define the dimensions of the nanowire
variable XEdge equal 100
variable YEdge equal 10
variable ZEdge equal 10
##################### Step 1: Initialization #####################
echo both
units metal
atom_style atomic
dimension 3
timestep ${dt}
boundary s s s
neighbor 2.0 bin
neigh_modify every 1 delay 10 check yes
################Step 2: Modeling and Setting ###################
Create a diamond lattice with <100> orientation
variable latconst_s equal 5.43
lattice diamond ${latconst_s} orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
Create a simulation box
region Simbox block 0 {XEdge} 0 {YEdge} 0 ${ZEdge} units box
create_box 1 Simbox
mass * 28.0855
Create a Nanowire
region Nanowire block 0 {XEdge} 0 {YEdge} 0 {ZEdge} units box
lattice diamond {latconst_s} orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
create_atoms 1 region Nanowire units box
##################### Step 3: Force Field #####################
pair_style tersoff
pair_coeff * * Si.tersoff Si
##################### compute eng #####################
compute eng all pe/atom
compute eatoms all reduce sum c_eng
##################### Create Regions #####################
variable Thick_LB equal 2*{latconst_s} #thickness of left boundary layer
region fixed_LB block INF {Thick_LB} INF INF INF INF units box
variable Thick_BS equal 2*{latconst_s} # Thicknes of bending stript variable lim_Bs equal {XEdge}-{Thick_BS} region Bending_strip block {lim_Bs} ${ZEdge} INF INF INF INF
group Gr_Bending_strip region Bending_strip
group Gr_fixed_LB region fixed_LB
group Gr_mobile subtract all Gr_fixed_LB
##################### Step 4: Minimization and Equlibriation #####################
min_style cg
minimize 0.0 1.0e-8 10000 10000
velocity all create ${T} 48284121 mom yes rot yes
fix 1 all nvt temp {T} {T} $(100.0*dt)
#fix 1 all nve
variable mypxxo equal “-pxx/10000”
variable mypyyo equal “-pyy/10000”
variable mypzzo equal “-pzz/10000”
variable mylxo equal “lx”
variable mylyo equal “ly”
variable mylzo equal “lz”
variable mytempo equal “temp”
fix def1 all print 10 "step {mypxxo} {mypyyo} {mypzzo} {mylxo} {mylyo} {mylzo} {mytempo}" file equil.{loadforce}.out screen no
dump 000 all custom 100 dumpequilibration.xyz id type x y z
thermo 100
thermo_style custom step temp pe ke etotal pxx pyy pzz lx ly lz
run ${iterequi}
unfix 1
unfix def1
undump 000
##################### Deformation #####################
reset_timestep 0
fix 1 Gr_fixed_LB setforce 0 0 0
fix 3 Gr_Bending_strip addforce 0 ${loadforce} 0
fix 4 Gr_mobile nve
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 tstep equal “step”
variable mypxx equal “-pxx/10000”
variable mypyy equal “-pyy/10000”
variable mypzz equal “-pzz/10000”
variable mylx equal “lx”
variable myly equal “ly”
variable mylz equal “lz”
variable mytemp equal “temp”
fix def1 Gr_Bending_strip print 10 “{tstep} {mypxx} {mypyy} {mypzz} {mylx} {myly} {mylz} {mytemp}” file displace.{loadforce}.out screen no fix def2 Gr_mobile print 10 "{tstep} {mypxx} {mypyy} {mypzz} {mylx} {myly} {mylz} {mytemp}" file stress.{loadforce}.out screen no
dump 1 all custom 100 dump.Si_bending_${loadforce} id type x y z
thermo 100
thermo_style custom step pxx pyy pzz lx ly lz temp etotal pe ke
run ${iterrun}
