Hello,
I am trying to find the Young Modulus values of the titanium crystal structure at 0 Kelvin. For stress calculation I use Clausius virial theorem or BDT (atomic) stress approach, both give the same result. In addition, since I work around 0 Kelvin in the stress calculation, the main effect comes from the potential part and there is no stress effect from the kinetic part, and when I include the temperature effect, there is no difference in the results.
When I use nvt or nve ensemble in the analysis, the results are in full agreement with the results in the literature, while when I use npt ensemble, there is a 35% decrease in the stress value in the direction I apply load and therefore in the young modulus value. Stress value appears to be zero because I reset the stress value in other directions.
I relax until the temperature and pressure are balanced. I take the barostat damping value as 1000xdelta(t) and the thermostat damping value as 100xdelta(t). If I increase the damping values, the system stabilizes more slowly and I have to run for longer. Thus, I reset the residual stresses, if any, on the part before starting the deformation. When calculating stress and strain in the analysis, I take the time average, the results do not change when I change the statistical averages.
Although I have tried the analyzes for different statistical situations and titanium crystal structures, I get correct results with the nve or nvt ensemble, and quite different results with the npt ensemble.
My intention was to simulate testing actually done in vacuum or atmospheric pressure, but for weeks now I can’t find the reason. I’ve included the script I used below, I’d be very grateful if you could help.
Kind regards,
Emel
Input file for uniaxial tensile loading of single crystal Titanium
#--------------- INITIALIZATION ---------------
units metal
dimension 3
boundary p p p
atom_style full
newton on
neighbor 0.1 bin
neigh_modify every 1 check yes
##--------------- ATOM DEFINITION ---------------
read_data Ti_alpha.data
##--------------- FORCE FIELDS ---------------
pair_style meam/spline
pair_coeff * * Ti.meam.spline Ti
#--------------- SETTINGS ---------------
timestep 0.001
variable ts equal 0.001
velocity all create 0.5 4928459 rot yes dist gaussian
#--------------- EQUILIBRATION & MINIMIZATION ---------------
dump 1 all atom 1 minimize.lammpstrj
minimize 1.0e-8 1.0e-8 10000 100000
undump 1
#-------------- RELAXATION ---------------
fix 1 all npt temp 0.5 0.5 0.1 iso 0 0 1 #or aniso
dump 1 all atom 500 relaxation_test_npt.lammpstrj
thermo 500
run 50000
undump 1
#--------------- COMPUTATION ---------------
compute 1 all stress/atom NULL
compute 2 all reduce sum c_1[1] c_1[2] c_1[3] c_1[4] c_1[5] c_1[6]
#--------------- VARIABLES ----------------
reset_timestep 0
variable Lx equal lx
variable tmp equal “lx”
variable L0 equal ${tmp}
variable Ly equal ly
variable Lz equal lz
variable Vol equal vol
variable up equal 2e-3
#--------------- DEFORMATION -----------------
unfix 1
reset_timestep 0
fix 1 all npt temp 0.5 0.5 0.1 y 0 0 1 z 0 0 1
fix 2 all ave/time 1 100 100 c_2[1] c_2[2] c_2[3]
fix per all ave/atom 1 100 100 c_1[1] c_1[2] c_1[3]
variable delta equal {up}*{tmp}
fix 3 all deform 1 x delta 0 ${delta} units box remap x
#--------------- Output strain and stress info to file for units metal
variable strain equal “(lx - v_L0)/v_L0”
variable STR equal “v_strain”
variable PXX equal f_2[1]/10000/v_Vol
variable PYY equal f_2[2]/10000/v_Vol
variable PZZ equal f_2[3]/10000/v_Vol
fix def1 all print 100 “{STR} {PXX} {PYY} {PZZ}” file strain_stress.def1.txt screen no
thermo 500
thermo_style custom step temp v_STR v_PXX v_PYY v_PZZ ke pe press lx ly lz vol
dump impactor all custom 500 dump.Ti.* id type x y z fx fy fz f_per[1] f_per[2] f_per[3]
dump 1 all atom 500 tensile_npt.lammpstrj
run 50000
SIMULATION DONE
print “All done”